Mathematical Behavior of MEFV Curves in Childhood Asthma and the Role of Curvature in Quantifying Flow Obstruction

Maximal expiratory flow-volume (MEFV) curves of pediatric patients are investigated using differentiation schemes and by computing their second derivative, d2V̇ /dV 2. Results show that spirometric tracings illustrate a characteristic well-defined behavior, where two distinct regions of the MEFV curve may be identified: (1) a concave profile during the initial expiratory maneuver, and (2) a convex profile over the greater lower region of the descending phase of the MEFV curve; this latter region is characterized by an approximately constant positive value of d2V̇ /dV 2 such that the descending MEFV limb may be captured by a quadratic function. Based on simple expiratory flow modeling, we show that d2V̇ /dV 2, and alternatively the local geometrical curvature κ(V), yield a measure of the relative degree of flow obstruction. In view of future clinical applications, we make use of an “average curvature index”, to assist clinician’s assessment of asthma severity, by quantifying curvature and summarizing global information in MEFV curves.


Introduction
Since their original introduction by Hyatt et al. [1], maximal expiratory flow-volume (MEFV) curves have become a widely used noninvasive lung function test, potentially sensitive to respiratory mechanics and flow limitation.In childhood asthma, spirometry is considered a standard tool for objective assessment of disease severity and regular reevaluation in the assessment of disease control [2,3].In particular, flow limitation is frequently quantified by the forced volume exhaled in the first second (FEV 1 ) to assess disease severity, and FEV 1 serves as the primary outcome in clinical studies [4].However, most asthmatic children have lung function values in the normal range independent of disease severity [4][5][6], and there exists a lack of correlation between FEV 1 measures and individual symptom scores in asthmatic children [7,8].Due to limitations in the effectiveness of FEV 1 , other lung function parameters which indicate airway obstruction have been calculated from the flow-volume (FV) curve including peak flow, forced expiratory flow at 25% of forced vital capacity (FVC), and forced expiratory flow at 75% of forced vital capacity [9][10][11].However, such parameters often lack repeatability and have not proved to be effective markers of flow obstruction [12].Indeed, classic lung function parameters, including FEV 1 , reflect single value points extracted from one location on the MEFV curve, and are then referenced to predicted values based on age, height, race, and gender.As a result, global curve information remains largely neglected, and interpretation of such indices, based on single point estimates, may be limited.
Given the limited usefulness of classic lung function parameters, some clinicians will commonly "eyeball" the MEFV curve to assess childhood asthma severity, where a concave pattern may be observed despite possibly normal FEV 1 values [10,13] and obstructive lung disease may be inferred upon the degree of curvature of the expiratory limb of the MEFV curve [12].Although widely spread, this approach remains however purely qualitative in nature and is based on a visual inspection of spirometric tracings.

ISRN Pulmonology
Alternative methods have been proposed to quantify the degree of curvature in FV curves and have attempted to encompass global curve information.Mead introduced a "slope-ratio" method based on the instantaneous tangent slope of the FV curve and sensitive to curvatures of the expiratory loop [14,15].A major drawback of his method remains the high sensitivity to local fluctuations in the FV curve.Other approaches have formulated a usable single index based on the shape of [16] or the area under [17] the MEFV curve but results remain largely sensitive to noise and neither method is implemented in clinical settings.More recently, a "curvature index" based on the mathematical definition of curvature has been introduced to quantify the curvilinearity of the FV curve [18].Although the method is promising, such "curvature index" reflects again a single extracted point from the FV curve, that is, the point of maximal curvature, and a somewhat arbitrary definition of the region of curvilinearty over the descending phase of the MEFV curve, that is, from 90% of peak expiratory flow (PEF) to 90% of FVC.Indeed, until present, the mathematical description of the descending phase of the MEFV curve, upon which flow limitation may be inferred, has relied largely on empirical approaches (e.g., best fitting equations), and the portion of the descending limb to be fitted has been arbitrarily selected with little mathematical or physiological rational [1,18,19].
In the present investigation, our aim lies in gaining further insight on the mathematical behavior of the descending phase of the MEFV curve which can be useful for quantifying flow limitation and subsequently for assessing childhood asthma severity.We first show that MEFV curves are best understood by implementing analytical differentiation schemes.In particular, the behavior of the second derivative of the flow with respect to expired volume, d 2 V/dV 2 , reveals two distinct regions of the FV curve in relation to its concavity/convexity and discerned by an inflection point (IP), where the condition d 2 V/dV 2 = 0 holds.The descending limb of the FV curve is predominantly convex (d 2 V/dV 2 > 0), and the function d 2 V/dV 2 exhibits a stable behavior at a nearly constant positive value for each patient.This latter result allows us to model the flow explicitly as a function of expired volume, namely, V = f (V ).Using a simple expiratory flow modeling approach, we qualitatively show that the second derivative, d 2 V/dV 2 , reflects a measure of airway resistance and hence flow limitation.Finally, we define an average curvature index (ACI), based on the concept of geometrical curvature, κ(V), which encompasses global curve information over the descending expiratory limb of the FV curve.We illustrate its potential clinical application for quantifying flow limitation in pediatric patients.The mathematical schemes presented here are tested for a range of MEFV curves obtained from pediatric patients with variable degrees of asthma severity.

Methods
The MEFV curves were obtained from pediatric patients who underwent pulmonary function tests using a Masterlab electronic spirometer (Jaeger GmbH; Würzburg, Germany).For each patient, we downloaded instantaneous flow-volume raw data (flow, V (V ), in liters/s and volume, V, in liters), characterized with uniform equidistant discrete points (measures were obtained with a uniform volume incremental step, ΔV, of 0.04 liter).Each patient curve was thus described by a discrete function, Vi = f (V i ), consisting of n uniformly spaced data points.A total of 10 patient curves were employed for the present study (4 females, mean (±SD) age, 10.1 ± 2.1 years).Their selection was initially based on a doctor's diagnosis, a previously measured reversible airway obstruction, and a positive skin prick test for at least one major inhaled allergen.
Noise, inherently present in the raw flow-volume data and presumably unrelated to mechanisms of flow limitation, limits the degree to which local features present in the MEFV curves may be characterized.These fluctuations may arise from instrument noise, rapid flow oscillations due to audible turbulence, or a possible lack of compliance from the patient.To obtain an estimate of a patient's MEFV curve free of noise, a prior smoothing of the raw flow-volume data is advantageous.Several methods have been suggested including ensemble averaging [20], which would require each subject to perform several spirometric maneuvers, or alternatively, direct digital filtering of the MEFV curve using a linear phase, finite impulse response, Hamming window filter [21].In the present study, we are specifically interested in obtaining continuous and smooth first, d V/dV, and second, d 2 V/dV 2 , derivatives of the MEFV curves.For such reasons, we implement a nonparametric smoothing spline algorithm [22] over the MEFV curve, which yields piecewise polynomials with smooth and continuous derivatives.The cubic spline estimator is constructed from minimizing the following function: where the smoothing parameter p is selected in the vicinity of 1/(1 + ΔV 3 /6) and the function, s(V i ), is the smoothed MEFV curve which consists of cubic polynomials between adjacent values, V i and V i+1 .Estimates for first, d V/dV, and second, d 2 V/dV 2 , derivatives of the MEFV curves can be based on the respective derivatives of s(V ), using numerical computation schemes.Based on Taylor series expansion, first, d V/dV, and second, d 2 V/dV 2 , derivatives are obtained from s(V ) using second-order finite difference schemes [23].
For a uniformly spaced function, first and second derivatives are then computed at any point, i, using central difference schemes: ( At the boundaries of the domain, first and second derivatives are similarly computed using right-and left-sided difference schemes, respectively, at the start and end point of the MEFV curve.

Results and Discussion
3.1.Behavior of MEFV Curves.Figure 1 plots the raw spirometric flow-volume data and the resulting smoothing spline approximations, s(V ), obtained for the 10 patients considered.The MEFV curves illustrate various degrees of curvature discernible by visual inspection of the descending limb.The descending limbs range from a quasilinear profile (Figure 1, second row, first column) to a strongly bent profile (Figure 1, fourth row, second column).Further examples include a patient curve which abruptly stops at the end of the expiratory maneuver (Figure 1, third row, second column) as well as a curve marked by a second expiratory peak (Figure 1, fourth row, first column), presumably due to patient noncompliance.We found that estimates of d V/dV and d 2 V/dV 2 obtained directly from the raw flow-volume data, Vi = f (V i ), using central difference schemes with no prior smoothing, were subject to large fluctuations which buried the baseline characteristic behavior of the second derivative of the MEFV curve.In contrast however, derivatives computed from the smoothing spline, s(V ), revealed the characteristic trend of d 2 V/dV 2 .Figure 2 illustrates corresponding second derivative, d 2 V/dV 2 , obtained from patients of Figure 1 and plotted as a function of expiratory volume V .As a general remark, all plots of the second derivative, d 2 V/dV 2 , follow a similar characteristic behavior, regardless of noise fluctuation or abrupt ending of the signal.That is, the function, d 2 V/dV 2 , may be described by two distinct regions.The first region, which begins at the start of the expiratory maneuver (V(t = 0)), corresponds to the monotonically increasing portion of the curve where d 2 V/dV 2 < 0. The interval over which d 2 V/dV 2 < 0 holds corresponds by definition [24] to the region over which the MEFV curve is concave.This region corresponds intuitively to the ascending limb of the curve which ends in the shape of a hat shortly beyond PEF.The second region corresponds to the interval over which the function, d 2 V/dV 2 , either remains steady at a constant positive value, d 2 V/dV 2 > 0, or converges towards a steady value after an overshoot or a period of rapid oscillation about the x-axis (see Figure 2).By definition, the region for which the condition d 2 V/dV 2 > 0 holds corresponds mathematically to the convex portion of the descending limb of the MEFV curve (similar to the shape of a bowl).In essence, each MEFV curve may therefore be mathematically summarized into two distinct regions: (i) a concave region (d 2 V/dV 2 < 0) described by the ascending phase of the FV curve and the "hat-shaped" region containing PEF, and (ii) a convex region (d 2 V/dV 2 > 0) defining the greater lower portion of the descending limb.
Following these results, we make two further assessments.Firstly, the concave and convex regions of the FV curve must be separated, by definition, at an inflection point (IP), mathematically defined by the volume, V IP , for which the condition (d 2 V/dV 2 ) IP = 0 holds [24].Following Hyatt et al. [1] and noting that the MEFV curve corresponds to the maximal expiratory flow that an individual can attain at any volume in his vital capacity (VC), we speculate that V IP may correspond to the degree of lung inflation at which the expiratory limbs of the isovolume pressure-flow (PF) curves begin to show maxima; however, this hypothesis would need further investigation beyond the scope of the present analysis.The second observation concerns the descending phase of the MEFV curve for which we observe a stable positive value of d 2 V/dV 2 .Integrating twice the second derivative, d 2 V/dV 2 , with respect to volume, V, over the interval ranging from V IP to the end point of the MEFV curve, V max , the descending limb may therefore be adequately described by the following quadratic function: where constants C 1 , C 2 , and C 3 must be solved for each patient.Moreover, from a visual inspection of Figure 2, we observe that for all patients, the second derivative stabilizes at a positive steady-state value in the vicinity of zero.Therefore, the absolute curvature of the convex region of the descending limb remains small; this will be shown shortly from the definition of the local curvature.In particular, for d 2 V/dV 2 −→ 0, (3) will reduce to a linear function, namely, V = C 2 V + C 3 (see example shown in Figure 7).
For each patient curve, we may then accurately describe the greater lower portion of the descending limb, by first identifying IP and then solving for the quadratic function (3) over the convex region, where d 2 V/dV 2 = 2C 1 (patient values of d 2 V/dV 2 = 2C 1 are indicated in subplots of Figure 2).In practice, the precise location of IP may be made difficult to identify for second derivative functions which exhibit oscillations about the x-axis (e.g., Figure 2, first row, second column).For such reason, the solution to (d 2 V/dV 2 ) IP = 0 is found by first running a lowpass digital filter over the function d 2 V/dV 2 [25].Beyond IP, the convex region is then fitted with the quadratic function (3) using a robust nonlinear least-square method with a Trust-Region algorithm [26].The identification of IP and resulting fitting schemes for the convex region of the descending limb are plotted in Figure 3.The quadratic model yields a good description of the convex region of the curves (mean (±SD) r 2 , 0.997±0.003;mean (±SD) RMSE, 0.042±0.003).Finally, we investigate whether the steady-state function, d 2 V/dV 2 , over the convex region of the MEFV curve (Figure 2) approaches values, d 2 V/dV 2 = 2C 1 , obtained from the quadratic fitting scheme of (3).To quantify the stability of the second derivative, d 2 V/dV 2 , over the convex region past IP, we compute at each discrete point, i, the average function, (d 2 V/dV 2 ) i , over the interval extending from V IP to V i : For brevity, results from (4) are plotted in Figure 4 for a restrained number of patients.Generally, values of (d 2 V/dV 2 ) i converge well over the convex region of the descending limb, towards d 2 V/dV 2 = 2C 1 .Such results further suggest that not only does the quadratic function of (3) describe accurately the descending phase of the MEFV curve past IP, but furthermore, the behavior of the convex region is indeed well captured by the nearly constant value of d 2 V/dV 2 = 2C 1 .

Interpretation of Second Derivative of MEFV Curve.
We attempt to qualitatively capture some of the physiological characteristics of the second derivative, d 2 V/dV 2 , and infer its potential usefulness to quantify flow obstruction.In the previous discussion, V denoted flow as observed in the MEFV tracings, namely, as a function of volume V .
Redefining both V and V as functions of time, t, we obtain flow, V (t), as a function of volume v = V (t), by rewriting flow as V (V −1 (v)), where V −1 is the inverse function of V .The function, V −1 , can be assumed to exist over an expiration phase because V is monotonically increasing in t.Making use of the formula for the derivation of an inverse function, , and applying the chain rule, we obtain where V (t) is the volume acceleration.Taking the derivative of ( 5) with respect to V and making use of differentiation rules for a quotient yields where ... V (t) is the third derivative of V with respect to time.Equation ( 6) therefore describes the second derivative, d 2 V/dV 2 , of the MEFV curve as a dynamic system, involving only derivatives of volume, V, with respect to time.For the sake of simplicity and to characterize qualitatively maximal expiratory flow dynamics solely over the descending limb of the MEFV curve, we restrain our analysis to a simple linear first-order, single-compartment lumped parameter model [27], analogous to an RC circuit, which may be modeled as where R is the airway resistance, C the lung compliance, and ΔP the change in transpulmonary pressure P(t).Following Abboud et al. [28], we assume for simplicity that the pressure, P(t), changes linearly with expiratory lung volume as and furthermore, the airway resistance, R(t), has a simple linear relationship to expired volume, namely, R(t) = R 0 + aV (t), where R o is the resistance at V (t) = 0 and a characterizes the positive slope of the linear function.By definition, lung compliance is given by C = dV in /dP, where V in = TLC − V, and TLC is the total lung capacity.Therefore, using (8), C may be approximated by a positive constant time-independent value.Using the model of ( 7), we derive expressions for both V (t), and ... V (t) (see Appendix).Substituting the resulting expressions for V (t), V (t), and ... V (t) into (6) yields upon simplification Based on (9), we can now qualitatively study the behavior of d 2 V/dV 2 using representative, yet physiological, parameters.For a 10-year-old patient [29], FVC is on the order of 2 liters, C ≈ 0.5 liter/kPa, and the transpulmonary pressure will vary approximately over an interval of 40 cm H 2 O such that parameter K ≈ 2 kPa/liter and is consistent with C = 1/K.For this generic pediatric patient, resistance is expected on the order of R o ≈ 0.7 kPa•s/liter [30].We vary the parameter a between 0.05 and 0.5 kPa•s, such that R(V) may be effectively increased for a severe obstructive case by dV 2 = 0.91 dV 2 = 0.93 dV 2 = 0.39 Corresponding plots of second derivative, d 2 V/dV 2 , versus expired volume, V, for the pediatric patients presented in Figure 1.
Coordinates are as shown in the top left corner subplot.Note that all plots illustrate a similar characteristic behavior: (1) an initial monotonically increasing profile where d 2 V/dV 2 < 0 which (2) settles towards a small finite positive steady-state value, with d 2 V/dV 2 > 0. Several examples illustrate an overshoot before converging towards a steady-state value (e.g.fifth row, second column), while others show a period of oscillation about the x-axis before reaching their steady-state value (e.g.first row, second column).Steady-state values of d 2 V/dV 2 are indicated in each subplot.a factor of up to 2.5 at V max .Figure 5 illustrates for such physiological parameters, the curve profiles for flow, V (t), and second derivative, d 2 V/dV 2 , when the airway resistance, R(V ), is increased through the variation of parameter a.Looking at Figure 5(a), we note that as the effective value of the resistance increases, the overall degree of curvature in the descending limb goes from a quasilinear profile (a = 0.05) to a concave shaped curve as a becomes large.For small values of parameter a, the resulting second derivative, d 2 V/dV 2 , remains approximately constant at a low positive value over the expired volume, with the curves for d 2 V/dV 2 running parallel to the abscissa (see Figure 5(b)).As resistance increases, curves of the second derivative are shifted towards higher values of d 2 V/dV 2 and gradually exhibit an apparent profile inversely proportional to the expired volume.Although our model remains an over simplification of real physical expiratory lung flow dynamics, our simple approach nevertheless illustrates qualitatively that a net increase in airway resistance yields a higher degree of curvature in the descending limb of the MEFV curve, and concurrently an increase in values of d 2 V/dV 2 .Moreover, (9) suggests that for relatively low degrees of flow obstruction, values of d 2 V/dV 2 may indeed remain approximately constant at low positive values, as seen for pediatric patients of Figure 2 over the convex region of the descending limb.Our present results therefore suggest that the relative degree of flow obstruction in childhood asthma may be measured with some confidence using the second derivative of the MEFV curve over the descending convex limb.

Curvature of the Descending Limb of MEFV Curves.
We have shown that the degree of curvature in the descending limb of the MEFV curve may be quantitatively assessed using values of the second derivative d 2 V/dV 2 .However, the direct geometrical interpretation of the dynamic system describing d 2 V/dV 2 in (6) may be somewhat nontrivial.Therefore, in an attempt to replicate the clinician's "eyeballing" of the MEFV curve [10,13], we introduce the concept of geometrical curvature and describe the local curvature, κ(V ), over the interval between V IP and V max , where the convex region is modeled by (3).Local curvature is given by the reciprocal of the radius of a circle fitted to the curve at a given point.The circle is then fitted such that it shares the same first and second derivatives with the curve at that point.This geometrical concept of curvature might describe more appropriately what meets the clinician's eye compared to the second derivative, d 2 V/dV 2 , of the MEFV curve.For a twodimensional curve written in the form V = f (V ), the local curvature, κ(V ), may then be written as follows [24]: From its mathematical definition, the local curvature, κ(V ), is a function of the second derivative, d 2 V/dV 2 , namely, κ(V ) ∝ d 2 V/dV 2 .Substituting into (10), first, d V/dV, and second derivatives, d 2 V/dV 2 , computed from the quadratic model (3), yield Figure 6 illustrates the resulting profile of κ(V ) obtained for various patient MEFV curves.For patients with low positive values of d 2 V/dV 2 , κ(V ) −→ 0 and remains approximately constant over the expired volume.As values of d 2 V/dV 2 increase, the κ(V ) curve is shifted gradually towards higher values further away from the abscissa and curves illustrate a rapidly increasing profile over the expired volume.Therefore, patients with relatively stronger flow obstruction (higher values of d 2 V/dV 2 ) will yield κ(V ) profiles with relatively higher values.Intuitively, higher values of κ(V ) correspond to stronger degrees of geometrical curvature visible by inspection of the MEFV curve.One may note that for spirometric tracings which exhibit a quasilinear profile in the descending phase, the second derivative converges towards d 2 V/dV 2 −→ 0, as illustrated in Figure 7. Keeping only firstorder terms in C 1 for MEFV curves with a quasilinear profile, (5) will then simplified to Intuitively, a perfectly linear profile in the descending limb of the MEFV curve has no curvature, as κ(V ) = 0 for C 1 = 0.
As we are interested in future clinical applications and the clinician's direct assessment of disease severity, we introduce a straightforward usable index, based on geometrical curvature, κ(V ), to quantify patient flow obstruction, compare relative levels of obstruction between patients, and obtain in a quantitative, reliable, and reproducible manner, the widely used qualitative "eyeballing" scheme.A mathematical translation of the "eyeballing" scheme will then correspond to the quantitative evaluation of an average curvature, κ, over the descending phase of the MEFV curve.In practice, the average curvature may then be defined as follows [31]: where the arc length is given by dl 1/2 dV and the length, L, of the convex descending phase of the MEFV curve is by definition In view of its potential clinical application [32], we refer to the average curvature, κ, as the "average curvature index" (ACI), which corresponds effectively to the integration of the local curvature, κ(V ), along the arc length, dl, of the descending limb, and averaged over the length, L, of the descending limb [32].Therefore, ACI encompasses global curve information over the descending phase of the MEFV curve and simultaneously delivers a simple usable index, comparable between patients and quantifying the degree of flow obstruction.Values of ACI are given for each patient subplot of Figure 3 from the computation of κ(V ) (12) followed by κ (13), using patient parameters, C 1 and C 2 , obtained from the quadratic fit (3).Further conceivable clinical applications of ACI may be foreseen, such as the quantitative assessment of the usefulness of inhaling bronchodilators in flow obstruction, as illustrated in Figure 8.As a final remark, given low values of d 2 V/dV 2 (see Figure 2) and following the definition of κ (13), typical values of ACI are expected to remain small, namely, ACI < 1.

Conclusions
Our results show that the behavior of MEFV curves is well understood using mathematical differentiation schemes to characterize maximal expiratory flow and, in particular, the shape of the descending phase of the FV curve.Furthermore, we have qualitatively shown that the relative degree of flow resistance is captured by the second derivative, d 2 V/dV 2 , of the MEFV curve.Therefore, flow limitation may be directly inferred by computing d 2 V/dV 2 over the convex region of the descending limb of the MEFV curve, or alternatively using the concept of geometrical curvature, κ(V ).From  a clinical perspective, the average curvature index (ACI) may reveal itself as a more practical and intuitive single index, in conjunction with classic lung function parameters (e.g., FEV 1 ), to assess asthma severity in pediatric patients.It is anticipated that algorithms computing ACI may be directly integrated within existing spirometric devices to obtain realtime assessment of asthma severity during lung function tests.

Figure 5 :
Figure 5: Physical model of maximal expiratory flow over the descending limb of the MEFV curve and resulting second derivative, d 2 V/dV 2 .Curves illustrate the influence of increasing the net resistance, R(V ), through parameter a, on (a) the descending phase of the MEFV curve and on (b) the second derivative of the MEFV curve.As the net resistance increases, the descending limb reaches higher degrees of curvature and the profile of d 2 V/dV 2 shifts from the abscissa towards higher values.Details on the choice of parameters are given in the text.

Figure 7 :Figure 8 :
Figure 7: Example of a patient with a distinct quasilinear profile in the descending phase of the MEFV curve.(a) Plot of second derivative, d 2 V/dV 2 , as a function of expired volume (b) MEFV curve with (--) resulting quadratic fit and (o) IP.Steady-state value of the second derivative, d 2 V/dV 2 , and ACI are, respectively, given in (a) and (b).