Efficient sampling for polynomial chaos-based uncertainty quantification and sensitivity analysis using weighted approximate Fekete points
EEfficient sampling for polynomial chaos-based uncertaintyquantification and sensitivity analysis using weightedapproximate Fekete points ∗ Kyle M. Burk , Akil Narayan , and Joseph A. Orr
Department of Biomedical Engineering, University of Utah, Utah, USA Department of Anesthesiology, University of Utah, Utah, USA Mathematics Department, University of Utah, Utah, USA Scientific Computing and Imaging Institute, University of Utah, Utah, USAAugust 12, 2020
Abstract
Performing uncertainty quantification (UQ) and sensitivity analysis (SA) is vital when devel-oping a patient-specific physiological model because it can quantify model output uncertainty andestimate the effect of each of the model’s input parameters on the mathematical model. By pro-viding this information, UQ and SA act as diagnostic tools to evaluate model fidelity and comparemodel characteristics with expert knowledge and real world observation. Computational efficiencyis an important part of UQ and SA methods and thus optimization is an active area of research.In this work, we investigate a new efficient sampling method for least-squares polynomial approxi-mation, weighted approximate Fekete points (WAFP). We analyze the performance of this methodby demonstrating its utility in stochastic analysis of a cardiovascular model that estimates changesin oxyhemoglobin saturation response. Polynomial chaos (PC) expansion using WAFP producedresults similar to the more standard Monte Carlo in quantifying uncertainty and identifying themost influential model inputs (including input interactions) when modeling oxyhemoglobin satura-tion, PC expansion using WAFP was far more efficient. These findings show the usefulness of usingWAFP based PC expansion to quantify uncertainty and analyze sensitivity of a oxyhemoglobindissociation response model. Applying these techniques could help analyze the fidelity of otherrelevant models in preparation for clinical application.
Performing uncertainty quantification (UQ) and sensitivity analysis (SA) is vital when develop-ing a patient-specific physiological model. UQ and SA can be used to quantify the relative effectof a model’s input variables (individually or interactively) on the variability of the response of aphysiological model.[57] By doing so, UQ and SA provide the basis for comparing model outputbehavior with expert knowledge and real world observations. Results from this comparison canthen be used to determine if the model is ready for clinical testing or if further adjustments arenecessary.Although UQ and SA were first used for assessing mathematical models of control and mechan-ical systems, they are also useful for assessing biological systems.[23] One method for ascertaininguncertainty and sensitivity with respect to varying input variables, when unable to compute the ∗ This work was supported in part by a Utah Space Grant Consortium Graduate Research Fellowship (KMB) and NIHaward U24-EB029012. a r X i v : . [ q - b i o . Q M ] A ug utput distribution directly, is a polynomial chaos (PC) expansion.[61] A PC expansion representsa model’s output based on values of random inputs.[16] This prediction is stochastic and, com-pared with simpler approaches such as computationally demanding Monte Carlo (MC) methods,can provide a more efficient representation of model output variability. A PC expansion can bemore efficient because after a one-time construction cost is invested, statistical quantities such asexpected values, variances, higher-order moments, and variable sensitivities can all be computedvia an inexpensive post-processing manipulation of the PC expansion.In PC expansion methods, the PC surrogate is formed by seeding the surrogate constructionalgorithm with model evaluations. A PC surrogate is a multivariate polynomial function of aspecified degree or order, with the model parameters serving as inputs to the polynomial. Withrespect to the number of degrees of freedom or expansion coefficients in the PC surrogate, it isrecommended that a multiplicative sampling factor f S be used so that the minimum number ofmodel evaluations is at least two times the number of expansion coefficients ( f S = 2).[19] Underthis requirement, the computational expense grows with the number of model runs which in turnincreases with the number of model parameters and the maximum polynomial order. Due to thislimitation, once the number of model parameters or maximum polynomial order surpasses a certainthreshold, PC expansion strategies may no longer be more efficient than MC methods.The computational expense for constructing PC surrogates can be reduced by either reducingthe number of model parameters, the polynomial order, or reducing the number of required modelruns. The number of model parameters can be reduced by first screening for important parametersbefore applying PC expansion.[19] This in effect reduces model dimensionality which in turn reducesthe required number of model runs, although it should be noted that screening may miss importantinteractions and understimate variability. A different approach, for reducing the number of requiredmodel runs, is to use a sparse polynomial chaos expansion.[6] This approach uses an iterativeprocedure that reduces the number of model runs by only adding basis terms that improve the PCexpansion.[5] This in turn reduces the number of required model runs by reducing the number ofunknown coefficients. One can investigate also decreases in the polynomial order of the PC model,but this effectively decreases the descriptive power of the surrogate, and a very small polynomialorder may blunt the surrogate’s expressivity and decrease the accuracy of quantities of interest. Thethird strategy of reducing the number of model runs has received much attention in the literature,and is usually concerned with developing model query strategies that minimize the requisite numberof model runs required to accurately construct a PC expansion. Our approach follows this thirdstrategy.Recently, Guo et al. have introduced a new efficient sampling method for least-squares poly-nomial approximation, Christoffel-weighted approximate Fekete points (WAFP).[24] Rather than amultiplicative sampling factor this sampling method may only require an additive factor. By usingan additive factor, this approach may decrease computational expense by directly reducing the re-quired number of model runs while still including all basis functions up to the maximum polynomialorder specified. The method is flexible with respect to the number of model evaluations generated,the polynomial order, the number of model parameters, and the probability distributions associatedto each parameter. In addition, the method requires only fairly simple linear algebra procedures tobe implemented.The purpose of this paper is to demonstrate the accuracy of a WAFP based PC expansion andits sensitivity indices when using an additive sampling factor. This is demonstrated by calculatingthe error of the PC expansion constructed using an additive sampling factor and by comparing itssensitivity indices with indices calculated using a multiplicative sampling factor as well as usingMC methods. After demonstrating the accuracy of a WAFP based PC expansion, this paperthen constructs and analyzes a PC expansion in order to analyze a model for describing changesin oxyhemoglobin dissociation response. The efficacy of the proposed approach on this examplemodel suggests that it can be used to efficiently perform stochastic analysis for UQ on more generalcardiovascular models.The remainder of this paper is organized as follows. The next section defines the notation thatwill be used throughout the paper. Section 3 explains the construction of a PC expansion andapproximation of its coefficients. In Section 4 the Christoffel-based weighted approximate Feketepoint approach is presented. Section 5 presents how uncertainty and sensitivity characteristicswere obtained from the PC expansion. Sections 6-7 detail the methods used for assessing a PCexpansion and its indices. Section 8 details the MC methods used for comparison while Section 9 escribes the oxyhemoglobin dissociation response model. The computational efficiency of WAFPbased PC expansion is demonstrated in Section 10 by applying the methods to a oxyhemoglobindissociation model. Lastly, Section 11 concludes the paper by discussing the WAFP procedure andits application to the dissociation model. Throughout this paper, an output of interest Y for a model f with k inputs in X = ( X , X , ..., X k )will be denoted as follows: Y = f ( X , X , ..., X k ) = f ( X ) , (1)where k is the number of independent stochastic input parameters X k . The output of interest Y isthe model’s stochastic output and is uncertain because the inputs in X are uncertain themselves.The stochastic inputs are contained by X . Here, we have defined the inputs, X , X , ..., X k , asuniform continuous variables which are mutually independent. The input parameter X i has proba-bility density ρ i . We assume that each input parameter X i has been linearly mapped to the domainΩ i = [ − , X i ∈ Ω i = [ − , Y is a physiological model, the domain Ω i and the distribution of the inputparameters will depend on the research question. When exploring the effect of measurement un-certainty on the variance of a model output, the domain should be defined using a measure ofvariation such as standard deviation or 95% confidence intervals. If known, the distribution of theinput parameters should match the known distribution of each type of measurement. The effect ofrealistic physiological ranges on a model output can also be explored. This analysis can help deter-mine the expected variance in a model output across the populations of both healthy and unhealthypatients. It can also determine the most influential input parameters on the model output. For thistype of analysis, the domain should be defined across the entire realistic physiological range and, ifknown, the input parameter distribution should match the distribution of experimental or clinicalmeasurements. If the distribution is unknown, the input parameters are usually assumed to have auniform distribution.Performing PC assumes that f ( X ) is smooth with respect to X across the sample space Ω X := × ki =1 Ω i . In other words, it assumes f ( X ) has many derivatives with respect to X .[22] Performingquality assessment of the PC expansion as detailed by Dubreuil et al. [21] and as implemented belowcan be used to verify the accuracy of the PC expansion approximation. If quality assessment testingshows that the PC expansion approximation is accurate than this assumption can be consideredvalid.Performing a quality assessment for PC surrogates can be used to determine when the approachin this article is useful and when it is not. If a PC emulator that accurately predicts the modelcan be constructed, then the procedures in this article will be beneficial for practitioners who wishto explore the model. However, some models are not smooth with respect to the input parameters X , and thus building a straightforward PC surrogate will not be effective. Alternatives to PCemulators are Gaussian process emulators[32], which may perform better in such situations. Sensitivity measures for the model were estimated using the Python software toolkit
Uncertain-SCI [1] which implements a WAFP-design-based PC method. In this section we provide detailsof the methods pyopoly uses for PC expansion construction and coefficient approximation. Later,Section 6 provides details of the methods pyopoly uses to quantify uncertainty (expected value,variance, and prediction interval) and analyze sensitivity (main and total sensitivity indices).
A truncated PC expansion of the output Y can be represented by the expression: Y ≈ (cid:101) f ( X ) = (cid:88) α ∈A c α Φ α ( X ) , (2) here c are the expansion coefficients, Φ are the polynomials, α are multi-indices with nonzeroentities, and A is the set of N p = (cid:0) k + pk (cid:1) multi-indices α for all polynomials with a maximal order p .We have constructed { Φ α } α ∈A to be L orthonormal polynomials with respect to density. Whenweighted by the joint density ρ x of the random variable X , i.e., they satisfy (cid:90) R k Φ α ( x )Φ β ( x ) ρ x ( x )d x = δ α,β , for all α, β ∈ A , where δ α,β is the Kronecker delta.For our model f , we further assume a uniform distribution for the input distribution ρ i of thecontinuous random inputs X i for simplicity. Having assumed this distribution, we use a particu-lar PC expansion basis, the Legendre polynomial family, for ease in computing and manipulatingsensitivities.[62] The construction of a PC expansion then proceeds by specifying a design of ex-periments over input variable space Ω x . The PC expansion’s maximal degree p of the total degreepolynomial approximation space is specified which results in a total of N p = (cid:0) k + pk (cid:1) polynomial basiselements. We note that our restriction to the uniform distribution in this paper does not limitthe applicability of our PC emulator approach: all of our procedures generalize to non-uniformdistributions. We construct the PC expansion with a weighted discrete least squares procedure. Least squaresapproximation typically requires a multiplicative sampling factor f MS to designate the number ofsampling points where N S = f MS N p . Least squares requires a minimum f MS = 1 sampling factor butgenerally f MS ≥ f AS where N S = N p + f AS . An additive sampling factor can be particularly efficient (i.e., requiring less data) when p (cid:29) k (cid:29) f AS of 10 is typically sufficient for the model example we consider. We note that our choice of f AS = 10is made for the types of models we consider. While we have experienced that the robustness of thisfactor is fairly uniform for the simulations we have investigated, such a oversampling factors shouldbe re-assessed for general situations.Explicitly, we assume that N S data points are available, (cid:110) x ( i ) , y ( i ) (cid:111) N S i =1 , where y ( i ) := Y ( x ( i ) ). Aweighted least squares procedure computes coefficients c α in (2) satisfying a standard least squaresproblem, { c α } α ∈A = arg min { d α } α ∈A ∈ R Np N S (cid:88) i =1 w ( i ) (cid:34) y ( i ) − (cid:88) α ∈A c α Φ α ( x ( i ) ) (cid:35) , (3)where the weights w ( i ) and the sample grid locations x ( i ) must be specified; the problem above is astandard least squares objective function, and can be solved explicitly via standard linear algebratechniques. The special distinction of the WAFP strategy is the choice of the weights and gridlocations. We describe this choice next. Once (3) is solved, the PC expansion surrogate f in (2) is determined and may be used forpredictive and analysis purposes. The method of WAFP determines a choice of the weights w ( i ) and sample locations x ( i ) in (3) [24]. We briefly describe this choice here. The weights are uniquelydefined by the sample locations, and are chosen as w ( i ) = (cid:34) (cid:88) α ∈A Φ α ( x ( i ) ) (cid:35) − . hus, once an experimental design for the locations x ( i ) is determined, then the weights are specified.The first step in determining the sample locations is construction of a large “training” grid, (cid:110) z (1) , . . . , z ( M ) (cid:111) ⊂ R k , M (cid:29) N S . (4)We will choose M ∼ KN p for some constant K , but at present we simply assume that this grid ispresent and discuss in Section 4.3 how this grid is generated.In what follows we detail how the WAFP design points are chosen. The procedure is completedin two steps: in the first step we choose N p sample points via one scheme, and in the second stepthe reamining N S − N p points are chosen via a different procedure. These steps are described inbrief in the following sections. The first step is identical to the original WAFP procedure [24]. Thesecond step is a modification of the original approach for oversampling, and is more computationallyefficient than the original approach. Although several other methods for experimental design exist,the WAFP strategy is simple and consistently produces stable PC expansion approximations inpractice[24]. N p points The true model Y is not evaluated on the training grid. Instead, the following model-independentmatrix is formed, B = ψ α (cid:16) z (1) (cid:17) ψ α (cid:16) z (2) (cid:17) · · · ψ α (cid:16) z ( M ) (cid:17) ψ α (cid:16) z (1) (cid:17) ψ α (cid:16) z (2) (cid:17) · · · ψ α (cid:16) z ( M ) (cid:17) ... ... . . . ... ψ α Np (cid:16) z (1) (cid:17) ψ α Np (cid:16) z (2) (cid:17) · · · ψ α Np (cid:16) z ( M ) (cid:17) ∈ R N p × M , (5)where ( α , . . . , α N p ) is any enumeration of the elements in A . The functions ψ α are defined as ψ α (cid:16) z ( i ) (cid:17) = √ w ( i ) Φ α (cid:16) z ( i ) (cid:17) . The matrix B is used to choose the sample locations { x ( i ) } N S i =1 . Precisely, a column-pivoted QR decomposition is performed on B , and the first N p column indices determined by this procedureidentify a size- N p subset of { z ( i ) } Mi =1 that defines the first N p sample locations. The selection isbased on analysis ensuring robust mathematical properties [9, 15]. Note that a column-pivoted QR decomposition simply performs column selection via greedy determinant maximization, so that inthe language of statistical design of experiments, this is a type of approximate D -optimal design.Note, however, that it is not a standard D -optimal design since a weighted determinant is maximzed. The remaining N S − N p samples are chosen via a slightly different procedure using a greedydeterminantal maximization, which is slightly different than the oversampling procedure describedin [24]. The previous section detailed the choosing of { x ( i ) } N p i =1 via a QR decomposition that iseffectively greedy determinant maximization. In this section we describe the choice of the set { x ( i ) } N S i = N p +1 , which is still greedy determinant maximization, but must be implemented slightlydifferently.The explicit iterative choice of x ( i ) for i > N p is given by x ( i ) = arg max z ∈{ z ( m ) } Mm =1 \{ x ( m ) } Npm =1 det( G ( i ) ( z )) , i = N p + 1 , . . . , N p + N S , (6)where G ( i ) ( z ) and B ( i ) are, respectively, i × i and N p × i matrices defined by G ( i ) ( z ) := B ( i ) ( z ) (cid:16) B ( i ) ( z ) (cid:17) T , B ( i ) ( z ) := ψ α (cid:16) x (1) (cid:17) · · · ψ α (cid:16) x ( i − (cid:17) ψ α ( z ) ψ α (cid:16) x (1) (cid:17) · · · ψ α (cid:16) x ( i − (cid:17) ψ α ( z )... . . . ... ... ψ α Np (cid:16) x (1) (cid:17) · · · ψ α Np (cid:16) x ( i − (cid:17) ψ α Np ( z ) ote that the explicit determinant above can be optimized more easily than the formulas suggest.If we define the N p -vector ψ ( z ) and the N p × N p matrix Ψ ( z ), ψ ( z ) := ψ α ( z ) ψ α ( z )... ψ α Np ( z ) , Ψ ( z ) := ψ ( z ) ψ ( z ) T , then we see that G ( i ) = Ψ ( z ) + G i , G i := i − (cid:88) r =1 Ψ ( x ( r ) ) , where G i does not depend on z . Thus, by the matrix determinant lemma,det( G ( i ) ( z ( j ) )) = det( G i ) (cid:20) ψ (cid:16) z ( j ) (cid:17) T G − i ψ (cid:16) z ( j ) (cid:17)(cid:21) . (7a)Therefore, maximizing the determinant in (6) is equivalent to maximizing the bracketed term above,which involves only a single matrix-vector and a single vector-vector multiplication, assuming G − i is stored from the previous iteration. In addition, once x ( i ) is computed from (6), then G − i +1 canbe computed from G − i without matrix inverses via the Sherman-Morrison formula: G − i +1 = G − i − G − i ψ (cid:16) x ( i ) (cid:17) ψ (cid:16) x ( i ) (cid:17) T G − i ψ ( x ( i ) ) T G − i ψ ( x ( i ) ) (7b)Therefore, the maximization (6) can be performed by simply using the update formulas (7) in aniterative fashion. This completes our description of how the WAFP points { x ( i ) } N S i =1 are chosen;our remaining task is to describe how the training grid is chosen. The quality of the WAFP nodes described in the previous subsections depends on the choiceof the training grid in (4). The WAFP procedures are determinant maximizing procedures thatuse linear algebraic operations. These operations in turn are most accurate when the matrix B iswell-conditioned. Thus, we seek to generate a good set of candidate samples so that the conditionnumber of B is as small as possible.In principle this can accomplished by a space-filling design, but this requires exploration of theentire parameter domain, which can be prohibitively expensive when k is large. To circumventthis, we leverage a least squares result regarding random sampling. Suppose that { z ( i ) } Mi =1 areindependent and identically distibuted random samples from the probability density ρ A ( x ) := 1 N p (cid:88) α ∈A Φ α ( x ) ρ x ( x ) (8)The density ρ A is called the induced distribution[39], and is a property only of the joint distribution ρ ( x ) and the polynomial index set A . With the z ( i ) distributed according to ρ A , then taking atraining set of size M satisfying for any C > B is well-conditioned with high prob-ability [15]. More precisely, the condition number inequality κ (cid:0) BB T (cid:1) ≤ − M − − C . Thus, an essential multiplicative behavior of M ∼ KN p for some constant K ensures that B is well-conditioned. This is the procedure we use to choose our training grid. Notethat efficient algorithms to generate random samples from a large class of induced distributions areavailable[39], but even a simple rejection sampling algorithm can easily generate samples from ρ A using samples of ρ x .For our WAFP training sets, we choose our candidate sets as random samples from ρ A butnote that alternative methods also generate good candidate sets[24]. Sampling points when usingthe WAFP procedure with maximal polynomial order p = 10 and sampling factor f AS = 0 for theOakley and O’Hagen[43] two-dimensional function are shown in Figure 1. .4 Alternatives to WAFP least squares In the sections above we have chosen the WAFP least squares strategy to compute the coefficients c α in the PC expansion (2). This is actually two choices: that of least squares approximation, and ofthe WAFP experimental design. We provide brief motivation of these choices in this section, alongwith discussion of alternative formulations that one can consider when constructing PC expansions.We do not give a comprehensive review of sampling strategies; more strategies and more completediscussions are available in the literature[42, 15].The least squares strategy, i.e., the formulation (3) defining the coefficients c α , could be re-placed by an interpolation strategy (if N S = N p ), or by more complicated methods such as sparseapproximation with compressive sampling [12, 46]. Our choice of the least squares formulation islargely motivated by the fact that simple, efficient, and stable algorithms can solve this problem,whereas more complicated algorithms can be required for alternative formulations. In addition,least squares formulations are intimately connected to statistical regression models, which have along history of theoretical and algorithmic development [20].Having selected least squares as the approximation technique, we can now discuss the particularsample design, i.e., the selection of the grid (cid:110) x ( i ) (cid:111) N S i =1 . In this article we cannot do justice to thevoluminous historical work on sampling and design of experiments for polynomial approximation;instead we focus on a discussion of recent advances and approaches, and seek mainly to appropriatelyalign and compare our approach with alternative existing methods. For least squares approxima-tions, we discuss below Monte Carlo methods, potential-theory-based least squares designs, anddesigns obtained via optimization. Perhaps the simplest approach for designing a sampling scheme is with randomization, whichproduces a random sample configuration. This approach is flexible, simple, and, generalizable.In this paper we have assumed that X is a random vector comprised of mutually independent,uniformly distributed random variables. Monte Carlo-type methods and algorithms apply to thiscase, but also apply straightforwardly without assumptions of independence or uniform distribu-tions. For random sample design PC methods, one constructs samples (cid:110) x ( i ) (cid:111) N S i =1 as independentand identically distributed (iid) realizations of the random variable X . One then computes thePC coefficients from (3) with w ( i ) = 1 for all i . The analysis in [14] reveals that, if N S is largeenough, then (3) is a well-conditioned regression design with high probability and computes a PCexpansion of near-optimal accuracy, i.e., computes a PC expansion whose L error is almost assmall as the best possible PC expansion, where “best” is measured in the L sense. The drawbackis that for most distributions of X , N S must be rather large compared to N P . For example, in ourcase of independent uniform random variables, we require N S (cid:38) N P using random, independentsamples. Note that this behavior of N S relative to N P guarantees convergence with high proba-bility, but often in practice one observes that N S = f MS N P for some constant f MS can empiricallyyield convergence.A slight variant of the random sampling algorithm, using importance sampling to bias thesampling distribution, can substantially improve convergence results. Such a weighted random-ized sampling procedure requires only that N S log( N S ) (cid:38) N P in order to obtain linear algebraicstability of the weighted least-squares problem (3) and convergence of ˜ f to the best polynomialapproximation [15]. The optimal biased distribution from which to choose samples is the induceddistribution shown in (8). The analysis then suggests that N S (cid:38) f MS N P samples are required (ig-noring logarthmic factors) to obtain accurate approximations [15]; in contrast the WAFP approachwe employ empirically requires only N S = N P + f AS samples [24]. However, we do not yet haverigorous theory that mathematically demonstrates the superiority of WAFP designs compared torandomized designs. We justify our choice of WAFP (in contrast to iid realizations), along withour recommended choice of f AS = 10 for the model considered in this manuscript in Figure 3. The second class of least squares designs that have become popular are those based on potential-theory concepts. (In the multivariate case these are pluripotential theory concepts.) In the frame- ork of least squares approximations, large polynomial degree asymptotics for the distribution ofan optimal sampling design are known[7]. More precisely, for a fixed polynomial degree k , one candefine a sampling design as optimal if it maximizes a particular weighted matrix determinant; thelarge- k limit of such designs must distribute according to a weighted pluripotential equilibrium mea-sure. Similar notions of such large-degree convergence to equilibrium measures are known in moreabstract settings [8]. This theory has inspired computational methods for least squares that exploitsuch a connection to pluripotential theory[40]. These pluripotential-theory-based approaches caneffectively generate good designs for large polynomial degree approximation if explicit forms forweighted pluripotential equilibrium measures is known. The disadvantages of these approaches arethat (i) the precise form of such weighted equilibrium measures are not explicitly known even inrelatively simple cases, making the development of computational methods difficult in general, and(ii) the optimality of such designs is true in the large degree limit; it is not clear how effective theseapproaches are for small-to-medium polynomial degrees. In contrast, the WAFP algorithm we ex-plore generates empirically good designs for an arbitrary polynomial subpsace, and does so withoutany knowledge of equilibrium measures. However, the theory for WAFP designs is still nascent andlacks some of the asymptotically-optimal properties of some other least squares designs. The last class of least squares designs we discuss are those generated by numerical optimization.Our WAFP approach fits best in this category. Such approaches attempt to compute a designby numerically optimizing a quality metric of the design. Such metrics can be moment matchingconditions [37, 33] or matrix determinants [53], but a wide variety of such metrics can be devised[41, 25]. More recently, compression techniques based on algorithms inspired by Carath´eodory’sTheorem and Tchakaloff’s theorem have been used to reduce a very large discrete set of samplesinto a smaller set of samples over which least squares procedures are stable [45, 34, 10]. Theseprocedures aim to generate parsimonious least squares designs that matches moments of a largediscrete sample set instead of a continous density. In contrast our WAFP approach matches themoments of the continuous density. In addition, the WAFP approach requires only standard linearalgebra procedures along with the explicit formulas in Section 4.3, whereas many of the previously-mentioned strategies require non-trivial optimization routines for implementation.
Uncertainty characteristics of expected value and variance are quantified directly using the PCexpansion. The prediction interval is estimated by approximating the probability density functionfor Y ( ρ Y ) using the PC expansion surrogate model and MC methods.The expected value and variance of ρ Y are obtained using the PC expansion. The expectedvalue is the PC expansion’s first expansion coefficient: E [ Y ] ≈ (cid:90) Ω x (cid:88) α ∈A c α Φ α ( x ) ρ x ( x ) d x = c . (9)The variance is the sum of squared expansion coefficients minus the first expansion coefficient: V [ Y ] ≈ E [( (cid:101) f ( X ) − E [ Y ]) ] = (cid:90) Ω x ( f ( x ) − c ) ρ x ( x ) d x = (cid:90) Ω x (cid:32) (cid:88) α ∈A c α Φ α ( x ) (cid:33) ρ x ( x ) d x − c . (10)Since percentiles and prediction intervals cannot be derived directly from the PC expansion of Y ,we apply the MC method to find them. A 95% prediction interval is created from the outputrealizations when evaluating the PC expansion for 10 , N s new input samples.We use Sobol indices to determine which parameters have little or no effect on the output Y andwhich parameters heavily influence Y .[56] Both main (first-order) and total indices are obtainedusing the PC expansion. Doing so requires approximating the total variance in Y by V [ Y ] ≈ (cid:88) α ∈A V [ c α Φ α ( X )] . (11) ain sensitivity ( S i ) is a unitless quantity between 0 and 1 that measures the relative variance acertain variable (without interaction) contributes to the total output variance V [ Y ]. For example, S measures only the variance that X imparts alone, and does not measure the variance contributedby, e.g., the pair ( X , X ) or any other pairs or higher-order tuples where X is involved. S couldrepresent the expected reduction in V [ Y ] that would be achieved if X were set to its standardphysiological value.[22] The main sensitivity index ( S i ) is calculated as: S i ≈ V [ Y ] (cid:88) α ∈A i V [ c α Φ α ] , where A i = { α | α i > ∧ α j = 0 ∀ j (cid:54) = i } . (12)Second-order sensitivity S ij is also a unitless quantity between 0 and 1 that measures the relativevariance that the interaction of two certain variables contributes to the total output variance V [ Y ].For example, S , measures the variance contributed by the pair ( X , X ). Higher order indicesalso exist, and these indices measure the variance contributed by higher-order tuples.[22]Total sensitivity ( S Ti ) measures the total relative variance a variable subset contributes, even ifother variables participate in this contribution. S T sums the first-order effect that X imparts alone( S ), along with the higher-order effects by, e.g., S , and any other pairs or higher-order tuplesinvolving X . If there is any variance imparted by the interaction between parameters, the sum of S Ti for all parameters will be greater than one. The total sensitivity index can be approximated by: S Ti ≈ V [ Y ] (cid:88) α ∈A Ti V [ c α Φ α ] , where A Ti = { α | α i > } . (13) S i is useful for ranking uncertain inputs and determining which inputs contribute to model outputvariability. S Ti is useful for analyzing the interaction between uncertain inputs and determiningwhich input interactions contribute to model output variability. If S i and S Ti − S i are near zero thisindicates that a particular input and its interactions do no contribute to model output variability.The input can then be set to a standard physiological value without affecting model fidelity. In order to properly evaluate WAFP-based PC expansion, the PC expansion in (2), as well asthe main and total sensitivity indices obtained from it, should be validated using a specific model(see Section 9). For this evaluation, an error measure for the PC expansion as well as the sensitivityindices should be calculated. This can be done by comparing the results for the PC expansion andsensitivity indices to reference values. For the analysis in this paper reference values were obtainedusing a WAFP-based PC expansion with a p = 10 maximal order PC expansion and f S,M = 2sampling factor. This assessment is performed in accordance with the previous work of Dubreuil etal. and Donders et al. which provide further detail in their manuscripts.[19, 21]
The error in the PC expansion (descriptive error) can be quantified using R , the coefficient ofdetermination (cid:15) R = 1 − R = (cid:80) N s i =1 (cid:16) y ( i ) − (cid:101) f P ( x ( i ) ) (cid:17) (cid:80) N s i =1 (cid:0) y ( i ) − y (cid:1) , (14)where y j is the sample mean of the realizations of Y j . The descriptive error (cid:15) R represents theresidual variance as a fraction of the total variance.The PC expansion can be validated using leave-one-out cross-validation, which computes thevalidation error using the same N s data points that are collected for construction of the PC expan-sion. The error in the PC expansion’s predictions (predictive error) can be quantified using Q , thevalidation coefficient Q = 1 − Q = (cid:80) N s i =1 (cid:16) y ( i ) − (cid:101) f ( − i ) P ( x ( i ) ) (cid:17) (cid:80) N s i =1 (cid:0) y ( i ) − y (cid:1) , (15)where (cid:101) f ( − i ) P represents a PC expansion constructed with all but the i th sample (i.e., with N s − N s PC expansion expansions are constructed to compute this metric. Thepredictive error (cid:15) Q represents the predicted variance as a fraction of the total variance. The relative error of the estimated (main and total) sensitivity indices S with respect to thereference sensitivity indices ˆ S can be used to quantify the accuracy of the sensitivity indices. Therelative error is (cid:15) δ = | S i − ˆ S i | ˆ S i , (16)where the index i = { , , ..., k } , because both main and total indices are considered. For theanalysis in this paper reference values ˆ S are obtained using a WAFP-based PC expansion with a p = 10 maximal order PC expansion and f S,M = 2 sampling factor.
To assess convergence, we set the multiplicative sampling factor f MS to 2 and the additive factor f AS to 10, both of which were WAFP PC expansions. Descriptive error ( (cid:15) R ) and predictive error( (cid:15) Q ) can then be calculated for 2 N PC expansions ( N multiplicative and N additive) generatedwith maximal PC expansion orders p = { , , , ..., N } . Descriptive error is calculated using (14)and predictive error using (15).PC expansions of a high maximal order p = N or less can be constructed to estimate theconvergence of error and variation in the sensitivity indices. Sensitivity indices are obtained fromeach of the 2 N PC expansions and are then evaluated for accuracy and precision compared to thereference analysis ( p = N + 1 and f MS = 2). Accuracy of the sensitivity indices is evaluated bycalculating relative error using (16). For UQ, the expected value, variance, and prediction interval are estimated from the evaluationsof Y using standard methods. For SA we used Saltelli’s alternative MC procedure as implemented inthe R software toolbox sensitivity [28]. This method used two sampling matrices of N samples eachdrawn independently of each other. For a four-dimensional model, k = 4 matrices are generatedfrom the two sampling matrices for a total of six matrices. Model evaluations are then calculatedfor each row in all six matrices. Further detail is provided in the literature for the interestedreader.[50, 22] We used Saltelli’s MC method with N = 100 ,
000 samples to quantify uncertainty and analyzesensitivity for comparison. To evaluate the convergence of Saltelli’s MC method, the relative errorof the sensitivity indices was also calculated using (16) for N = { , , , } . Referencevalues were obtained using a WAFP-based PC expansion with a p = 10 maximal order PC expansionand f S,M = 2 sampling factor. Oxyhemoglobin dissociation model
The oxyhemoglobin dissociation curve (ODC) describes the relationship between the partial pres-sure of oxygen in the blood (P o ) and the percent of hemoglobin saturated with oxygen (SHbO ).This relationship varies from patient to patient.[26] The position not the shape of the ODC varies,where for a given P o hemoglobin increases oxygen uptake as the curve shifts to the left and de-creases uptake as the curve shifts to the right.[26, 51, 31] P , the P o at which hemoglobin is 50%saturated (Figure 2), is the inflection point of the ODC and thus is commonly used to representthe position of the entire ODC. Determining a patient-specific P is helpful as it can help predictthe SHbO at which the patient’s oxyhemoglobin saturation will transition from declining slowlyto rapidly.P is determined by four factors (three chemical and one physical).[60] These four factors arethe concentration of the hydrogen ion (pH), temperature ( T ), the partial pressure of carbon dioxide(P co ), and the concentration of 2,3-diphosphoglycerate (DPG).[44] These factors influence P either directly or indirectly through each other. Many groups have studied these factors to createmodels capturing each parameter’s influence on P .[31, 60, 52, 55, 11, 35] These groups haveassessed the effect of these parameters as well as their synergistic effects using experimental dataand computational methods in order to model transport and exchange of oxygen in the lungs andtissue.[17] Most recently, Dash, Korman, and Bassingthwaighte described and evaluated an equationof state which incorporates the influence of all four variables (both individually and in combination)on P .[18]A UQ and SA of Dash’s equation of state, which has 4 input parameters and 4 polynomialequations, could formally quantify the uncertainty of P and its sensitivity to variations in pH, T ,P co , and 2,3-DPG, thus providing a means for ranking these variables and their interactions frommost to least significant.One utility of UQ and SA is that they can be used to analyze the fidelity of a model by comparingmodel characteristics to expert knowledge and real world observation. Using UQ and SA for thispurpose can help determine if a model is prepared for clinical evaluation. If the model is not ready,UQ and SA can be used to determine which characteristics of the model’s output do not align withtheory and can provide insight and guidance for aligning model characteristics with theory.The oxyhemoglobin dissociation model is suitable for highlighting this particular utility of UQand SA. Oxyhemoglobin dissociation was first described in the early 1900s and has been discussedextensively in the literature. The primary physical and chemical factors effecting oxyhemoglobindissociation are well-established. The influence of these factors on oxyhemoglobin dissociation hasbeen studied both experimentally in the laboratory as well through observing pathophysiologyof adjustments to different climates, situations, and diseases. This broad scope of knowledge onoxyhemoglobin dissociation makes the oxyhemoglobin dissociation model an ideal candidate fordemonstrating the benefits of applying UQ and SA to evaluate model fidelity in preparation fortransitioning to clinical use. Dash et al. analyzed Buerk and Bridges’ model and experimental data to determine curve-fitpolynomials for characterizing the relationship between P and changes in pH, T , P co , and 2,3-DPG [concentration].[18, 17] They determined these equations by varying one variable while theother three were fixed at standard values. Their results showed that the best-fit polynomial for therelationship between P and changes in each variable is given by: P , ∆ Q i = P S + (cid:88) j =1 Q ij ( Q i − Q S ) j , (17)where the vector Q = { pH, T, P CO , [ DP G ] } contains the physiological levels of each of the fourvariables. Q ij represents the coefficients for each of these four variables. These coefficients aregiven in Table 1. The quantities with superscript “S” denote the standard physiological value or that variable (see Table 2). Previous studies [11, 18, 31, 52, 55] have shown when combiningsimultaneous changes in more than one variable, P is given by: P = P S
50 4 (cid:89) i =1 (cid:16) P , ∆ Q i P S (cid:17) (18)Combining (17) and (18) gives the overall description for how P changes in response to varyingphysiological conditions: P = P S
50 4 (cid:89) i =1 (cid:16) P S + (cid:80) j =1 Q ij ( Q i − Q S ) j P S (cid:17) . (19)P is the output of interest for the oxyhemoglobin dissociation model. It is uncertain due tointerpatient differences in pH, T , P co , and 2,3-DPG which are uncertain themselves due to mea-surement uncertainty and stochastic variation over time and space within individuals. The analysisfor this portion of the experiment focuses on the uncertainty and variability across the entire pop-ulation. When subjected to uncertainty, we can express the four-parameter mathematical model( f P ) described by (19) in black-box form as: P = f P ( pH, T, P CO , [ DP G ]) = f P ( X ) , (20)where P is the stochastic output representing P which in this case is a scalar value. P isthe random field of interest and thus is a function of X which we have denoted by f P ( X ). Thestochastic inputs are contained by X = [ pH, T, P CO , [ DP G ]].We have defined the physiological sample space (Ω i ) for each input in Table 3. These are therealistic physiological ranges for each parameter. By analyzing these ranges, we have considered thecase where each parameter’s uncertainty is across the entire realistic physiological range. Therefore,the physiological range was the 100% confidence interval for that given input. This in effect evaluatesthe uncertainty of the output due to the inherent physiological variability of each variable in bothhealthy and unhealthy patient populations. We analyzed the role of pH, T , P co , and 2,3-DPG on the uncertainty of P by constructinga WAFP based PC expansion. (19) provides a model which maps the four input variables pH,P co , T , and 2,3-DPG [concentration] to the value of P . Constructing a PC expansion providesa surrogate model from which uncertainty and sensitivity can be extracted directly.We defined the inputs, pH, T, P CO , [ DP G ], as uniform continuous variables which are mutuallyindependent. We used an additive sampling factor f AS = 10 and maximal polynomial order p = 3which we confirmed were acceptable with the pre-processing analysis performed in Sections 6-7 andin Appendix A. This resulted in a PC expansion of N p = (cid:0) (cid:1) = 35 basis elements constructedfrom N s = 35 + 10 = 45 model runs. The model’s 95% prediction intervals were then estimated byevaluating the PC expansion using 10 , N s = 450 ,
000 new input samples of each input parameter.For comparison, we calculated model output uncertainty and sensitivity indices using MC methodsas detailed in Section 8.
We evaluated the output of the oxyhemoglobin dissociation model to evaluate the utility of usingPC expansion for this model, as opposed to MC methods. We varied the value for each parameter,one at a time, within the physiologically realistic parameter space. We did this both while holdingthe other input parameters to their standard values and while allowing the other input parametersto vary randomly within their sample space. able 1: Value of coefficients Q j in the P model (Eq. 17) from Dash et al.[18] j pH j P CO j [ DP G ] j T j − . · . · − . · . · . · . · − . · . · − − . · . · . · . · − Table 2:
Standard values of the oxyhemoglobin dissociation model
Parameter Description Value UnitpH S Standard pH in plasma 7.4 Unitless T S Standard temperature of blood 37 o CP co ,S Standard partial pressure of CO in blood 40 mmHgDPG S Standard 2,3-DPG concentration in RBCs 4.65e -3 MP ,S P at pH S , T S , P co S , and DPG S Table 3:
Uncertainty ranges (Ω i ) tested for the oxyhemoglobin dissociation model Nr. Symbol Uncertainty Range Unit1 pH 6.0-8.2 Unitless2 T 20-43 o C3 P co -3 -10e -3 M Table 4:
Results of the UQ of P for MC and PC with uncertainty ranges across the entire physiological range. Method Mean Standard deviation 95% Prediction interval[mmHg] [mmHg] [mmHg]MC 33.92 23.35 [7.33, 94.73]PC 34.26 23.15 [7.58, 94.16]
Table 5:
Polynomial model sensitivities
Sensitivity Analytical Order of PC expansionindex value p = 3 p = 4 p = 5 p = 6 S S S S S S S N p
20 35 56 84 N s
30 45 66 94 able 6: Ishigami function sensitivities
Index Analytical Order of PC expansionvalue p = 3 p = 5 p = 7 p = 9 S S S S S S S S T S T S T N p
20 56 120 220 N s
30 66 130 230
Table 7:
Full G-function sensitivities
Index Analytical PC-basedsolution solution ( p = 2) S S S S S S S S S T S T S T S T S T S T S T S T Table 8:
Reduced G-function sensitivities
Index Full model Reduced Model(eight parameters)Analytical p = 3 p = 5 p = 7 S S S S S S S S S S S T S T S T S T N p
35 126 330 N s
45 136 340 After constructing PC expansions for p = { , , , ..., } , we assessed the quality of the PCexpansions using a validation test set. The validation test set is generated independently of theexperimental design and WAFP training set. We recall that from Section 4.3 that the WAFPtraining set is generated from the density ρ A in (8). The validation set is generated as independentrandom samples from the density ρ . The discrepancy on this validation set was computed for thePC expansion prediction of P compared with the true value.First we assess both our choice to perform sampling via the WAFP approach, and our recom-mendation of choosing N S = N p + f AS samples, with f AS = 10. Figure 3 shows the behavior of WAFPsampling both as f AS is varied for N S = N p + f AS samples, and as f MS is varied for N S = f MS N p samples. The results demonstrate that the error depends mildly on f AS once f AS (cid:38)
5. For safety,we therefore choose f AS = 10. If we instead consider a sample count N S = f MS N p , the resultsdemonstrate that one can also attain attractive error results for say f MS (cid:38) .
3, but this generallyresults in many more samples, so that the additive recommendation N S = N p + 10 achieves similarerror with substantially reduced cost.Figure 3 also compares the WAFP strategy against error for an alternative PC constructionbuilt from N S iid Monte Carlo samples. The results show that the W AF P strategy consistentlyoutperforms the Monte Carlo strategy.Predictive ( (cid:15) Q ) and descriptive errors ( (cid:15) R ) for the PC expansion are shown in Figure 4. Amaximum polynomial order of at least p = 3 was required to obtain errors less than 10 − . PCexpansions with order p ≥ (cid:15) Q (cid:28)
1) sampling power for all q. Discrep-ancy between the PC expansion description and the true model was low ( (cid:15) R (cid:28)
1) with maximalpolynomial order p ≥ Relative error ( (cid:15) δ ) of the sensitivity indices for the additive and multiplicative PC expansionsare shown in Figures 5 and 6. In general, both PC expansions required a maximal polynomialorders p = 3 to obtain a sensitivity error less than 10 − . Although MC analysis performed ordersof magnitude more model runs, it did not obtain a sensitivity error less than 10 − . The estimated distribution of P was similar for both methods (Figure 7). The mean, standarddeviation, and 95% prediction intervals for MC and PC methods are shown in Table 4. Whenmodel inputs are sampled uniformly across a realistic physiological sample space, P does not havea normal distribution. The distribution is instead asymmetric with positive skewness. The sensitivity indices calculated using PC and MC methods correspond well (Figure 8). Theparameter ranking for influence on model output was the same when using both methods. Forthe concentration of hydrogen ions (pH) and temperature ( T ), both methods also reported totalsensitivities that were different from main sensitivities, indicating a substantial interaction betweenthe input variables.The concentration of hydrogen ions (pH) was estimated as the most influential variable, followedby temperature ( T ), the partial pressure of carbon dioxide (P co ), and the concentration of 2,3-DPG ([DPG]). For these estimates, the concentration of hydrogen ions ( S T ≈ .
72) had more thantwice the influence on total output variance than the second-ranked variable, temperature ( S T ≈ . > .
99, accounting for nearly all of the variance in the model’s output. Mutual interactions among aired variables accounted for 8 .
64% of variance in P , 7.8% of which was attributed to pairwiseinteractions between pH, T , and P co . One of the interactive effects had more influence on un-certainty of P than individual variables. The interaction between the pair (pH, T ) was rankedhigher ( S ≈ .
07) than both P co ( S ≈ .
03) and 2,3-DPG ( S ≈ . Although for variables the extreme P values occurred at the ends of each range, the dependencyof the equation of state output is non-linear but smooth across the entire input parameter space(Figure 9). The smoothness of the map therefore makes a PC expansion-based global UQ andSA an appropriate analysis for this model.[61] The smooth shape of the dependency justifies theuse of a global PC expansion as opposed to MC methods. Since the relationship is smooth, MCsimulation techniques would not be as effective or accurate in providing a sample representativeof the parameter space. On the other hand, a WAFP based PC expansion efficiently provides anaccurate emulator of smooth models. −0.010 −0.005 0.000 0.005 0.010x −0.010−0.0050.0000.0050.010 x WAFP 2-D sampling grid
Figure 1:
Two-dimensional grid generated using the weighted approximate Fekete points procedure.
11 Discussion
In this paper, we have proposed a procedure for computing the output uncertainty and sen-sitivity indices of physiological models with reduced computational cost. This procedure reducescomputational cost by constructing an inexpensive surrogate model using PC expansion. This pa-per focuses on reducing the computational cost of constructing a PC expansion by proposing adeterministic least squares regression sampling procedure, Christoffel-based weighted approximateFekete points (WAFP).The results show that using an additive sampling factor performs comparably to a multiplicativesampling factor. For both additive and multiplicative methods, a maximal polynomial order ≥ − indicating that the order of polynomialnecessary was independent of sampling strategy.. However, an additive sampling factor is moreefficient since it requires fewer model runs and can achieve the same accuracy using a fewer numberof points. This becomes particularly important when p (cid:29) k (cid:29)
1, where p is the PC expansion’smaximal degree and k is the number of independent stochastic input parameters. igure 2: Changes in oxyhemoglobin dissociation curve shape and position. The black solid line indicates how thepartial pressure at which 50% of hemoglobin are saturated with oxygen (P ) changes. P , and thus the oxyhemoglobindissociation curve, varies from patient to patient where oxygen affinity and P have an inverse relationship (oxygeninfinity increases as P decreases and vice versa). N S f AS L v a li d a t i o n e rr o r PC medianMC medianPC 5/95 quantile regionMC 5/95 quantile region 35 40 45 50 55 60 65 70Sample count N S f MS L v a li d a t i o n e rr o r PC medianMC medianPC 5/95 quantile regionMC 5/95 quantile region
Figure 3: L validation error on the P model for the WAFP PCE procedure proposed in this manuscript (“PC”,red) versus a PC construction that uses iid Monte Carlo samples (“MC”, blue). Left: N S = N p + f AS samples. Right: N S = f MS N p samples. The validation error is computed over 5000 iid samples generated independently of the samplingprocedures. Since the experimental design is random for both of these sampling procedures, quantile regions were com-puted using 100 realizations of each kind of PCE construction. The PC procedure considered in this manuscript performsnotably better than a standard MC approach. ncreasing maximal polynomial order1 3 5 7 9 11 Figure 4:
Predictive ( (cid:15) Q ) and descriptive ( (cid:15) R ) errors for the PC expansion of the oxyhemoglobin dissociation model. −10 −8 −6 −4 −2 S e n s i t i v i t i e s e rr o r [ - ] pH PCE f MS PCE f AS MC T Number of runs [-]10 −10 −8 −6 −4 −2 S e n s i t i v i t i e s e rr o r [ - ] CO Number of runs [-]
DPGMain sensitivity (S i ) error for the P50 model Figure 5:
Relative error ( (cid:15) δ ) in the main sensitivity ( S i ) indices for P obtained using PC expansions with additive andmultiplicative sampling factors and with increasing polynomial orders. The same relative error is plotted of the referenceMC analysis for N = { , , , , } samples per parameter. −10 −8 −6 −4 −2 S e n s i t i v i t i e s e rr o r [ - ] pH PCE f MS PCE f AS MC T Number of runs [-]10 −10 −8 −6 −4 −2 S e n s i t i v i t i e s e rr o r [ - ] CO Number of runs [-]
DPGTotal sensitivity (S Ti ) error for the P50 model Figure 6:
Relative error ( (cid:15) δ ) in the total sensitivity ( S Ti ) indices for P obtained using PC expansions with additive andmultiplicative sampling factors and with increasing polynomial orders. The same relative error is plotted of the referenceMC analysis for N = { , , , , } samples per parameter. Implementing a WAFP based PC expansion can yield results that are the same as MC methods.However, the computational cost is much lower with WAFP based PC expansion. One reason thecomputational cost is lower is because WAFP based PC expansion is able to use an additive samplingfactor. Using an additive sampling factor means that WAFP based PC expansion sampling sizegrows less as the maximal polynomial order and model dimensionality grow.In this paper, we have used uncertainty quantification to analyze a four parameter model fordescribing changes in oxyhemoglobin dissociation as represented by P . This uncertainty analysishas included analysis of the forward uncertainty propagation of P and of the sensitivity of P to input parameters pH, T , P co , and 2,3-DPG as expected using the oxyhemoglobin dissociationmodel. This quantification has shown that the expected distribution of P when using the modelis positively skewed with the 95% prediction interval for P ranging from 8 to 94 mm Hg. Thisanalysis has shown that, as implemented in the model, P is most sensitive to perturbations ofpH, sensitive to T , and insensitive to P co and 2,3-DPG. Therefore, only pH and T are interestingfor modeling the uncertainty in P due to physiological variability. Based on our simulationresults using the WAFP procedure, modeling P co and 2,3-DPG does not contribute significantlyto uncertainty in P .The highest ranked parameter for sensitivity of Dash’s equation of state aligns with theoreticaland expert knowledge. pH’s large influence on oxyhemoglobin dissociation, called the Bohr Effect,plays a major physiological role in transporting oxygen. The Bohr Effect is the mutual interactionof oxygen binding and hydrogen ion binding of oxyhemoglobin [54]. This mutual interaction allowsa decrease in pH near tissue capillaries to increase delivery of O at the tissue level [18, 48]. Thusthe high influence of pH on P plays a critical role in hemoglobin’s ability to unload O to tissues.The pairwise interactions between influential parameters shown using sensitivity analysis shouldreflect theoretical and expert knowledge on how these parameters influence each other. Dash’s equa-tion of state showed substantial interaction between pH and T . This result aligns with experimentalresults from studies conducted separately from the studies used to establish the equation of state.The interaction between these parameters has been described previously and shown to have a linear igure 7: Numerical evaluation of the distribution for P using Monte Carlo (MC) and polynomial chaos (PC) tech-niques. MC with N s = 100 ,
000 is shown on top and PC with order p = 3 and f AS = 10 is shown on bottom. The centervertical thick line denotes the estimated value µ ( P ). The dashed lines to either side denote µ ( P ) ± σ ( P ). Thedot-dashed lines represent the 95% prediction interval. S en s i t i v i t y [ - ] S i : main indexS iT : total index PC MC CO [DPG]pH T PC MC PC MC PC MC
Figure 8:
Main ( S i ) and total ( S Ti ) sensitivity Sobol indices for pH, T , P co , and DPG with P as the output ofinterest. Estimates from polynomial chaos (PC) method are shown with the left two bars of each group. Estimates fromthe Monte Carlo (MC) method are shown with the right two bars. For each method, S i is shown on the left and S Ti onthe right.
50 100 150 P (mmHg) I npu t p H [ - ] pH P (mmHg) I npu t T e m pe r a t u r e [ o C ] T P (mmHg) I npu t P C O [ mm H g ] CO P (mmHg) I npu t [ D P G ] [ M ] -3 [DPG] Figure 9:
Graphical portrayal of the non-linear output behavior for the oxyhemoglobin dissociation model which hasfour interacting model inputs, pH, T , P co , and 2,3-DPG. The solid black line indicates the value of each individualparameter across its physiological range with the three other parameters at standard values. The gray dots represent theoutput P = f ( pH, T, P CO , [ DP G ]), plotted as the individual parameter against the output, with the values of theother parameters varying across their physiological range. As shown, the response in-between each parameter’s endpointsis non-linear, indicating that the use of a PC expansion would be advantageous as opposed to Monte Carlo methods relationship [49].Dash et al. have developed an accurate equation of state for determining SHbO under variousphysiological conditions.[18] Dash’s equation of state has been well-validated on a large set of ex-perimental data from separate experiments available in the literature.[29, 38, 2, 27, 36, 47] Dash hasused these data to characterize the dependence of P on pH, T , P co , and 2-3-diphosphoglycerate(2,3-DPG). This equation of state is the most extensively validated using independent data andthe most accurate at the time of this writing. When implementing their equation of state, theirmodel provided accuracy greater than 99.5% over the whole saturation range when compared withexperimental data.Despite these efforts and this accuracy, the oxyhemoglobin dissociation model does not entirelyalign with theoretical and expert knowledge. The main ( S ≈ .
01) and total ( S T ≈ .
02) sensi-tivities for 2,3-DPG do not support well-established knowledge that 2,3-DPG alters P both bybinding to hemoglobin and by altering pH [58, 3, 13, 59]. Analysis of the oxyhemoglobin dissociationmodel shows that for the model the influence of 2,3-DPG is negligible compared to the influence ofpH, T , and P co . One explanation for the lack of sensitivity to 2,3-DPG could be that variationsin the normal physiological range do not change the output even though the unphysiological rangeclose to 0 would drastically affect the oxygen release from hemoglobin. However, without 2,3-DPGalmost no oxygen would be released from hemoglobin. This implies that 2,3-DPG, either directlyor indirectly, does in fact influence oxyhemoglobin dissociation.Despite conservatively testing a broad range of P co , Dash’s equation of state shows minimalinteraction between pH and P co ( S pH,CO ≈ . co inhuman blood is well-established and described by the Henderson-Hasselbach equation [30]: pH = pK a,H CO + log (cid:16) [ HCO − ][ H CO ] (cid:17) . (21)This interaction is an essential part of the bicarbonate buffer system for maintaining pH in the lood. Without it the body would be unable to maintain acid-base balance.Up to this point, the focus when developing oxyhemoglobin saturation models has been on theaccuracy of the model across the full range of P o . The analysis in this paper shows the importanceof also performing UQ and SA to analyze the fidelity of a model. Despite high accuracy whenanalyzed using experimental data, the Dash model only partially aligns with clinical observation andexpert knowledge on changes in oxyhemoglobin saturation. Performing UQ and SA has highlightedthe characteristics of the model which require adjustment before model fidelity can be obtained.Two of these areas include the influence on P of DPG individually and of the interaction betweenpH and CO . Only after model characteristics are improved in these two areas will the Dash modelbe able to properly simulate changes in P due to different levels of these effectors.In their paper, Dash and Bassingthwaighte initially concluded that pH and T have the mostinfluence on O binding.[17] After refining their model they concluded that pH significantly influ-ences O binding compared to the smaller effects of P co , 2,3-DPG, and T .[18] UQ and SA was notincorporated into Dash et al.’s studies and thus they did not reach their conclusions using formalquantitative computational methods nor did they provide a description of the interactions amongvariables. We have shown the value of UQ and SA by providing a precise quantitative descriptionof the uncertainty of P and a ranking of the sensitivity of P to all four variables. We have alsoshown the value of UQ and SA by providing a formal analysis of the interactions among variables.The Dash model is interesting to analyze because it has two influential (pH, T ) and two nonin-fluential parameters (P co , DPG). It is also interesting because two of the variables interact. Theseinteresting aspects make the Dash model a good candidate for future benchmarking of uncertaintyquantification techniques.Many other models exist which describe oxyhemoglobin dissociation. In future research, Buerkand Bridges’ model [11], which Dash and Bassingthwaighte based their initial model on, could beevaluated for comparison. Further, the methods described in this paper could also be performedon models which incorporate different subsets of the four model input parameters analyzed here.Analysis of these models could then be compared with the analysis in this paper to determine whichmodels align best with real world observation.We have assumed a uniform distribution of pH, T , P co , and 2,3-DPG across the parameterspace. Right now there are no results providing evidence of a different distribution. With a uniformdistribution, any value of pH, T , P co , and 2,3-DPG is equally likely to be sampled. Should datadescribing the distribution of one or all of these parameters become available in the future, thatdistribution could be used in a PC expansion, and an analysis identical to the one presented herecould be used to study uncertainty and sensitivity, with appropriate modifications to the inputdistribution.The focus of this paper has been on implementing a WAFP based PC expansion procedure forreducing computational cost. Other methods for reducing computational cost have been describedpreviously including input parameter screening prior to constructing a PC expansion as well assparse PC expansion.[19, 5] Although this paper has focused on WAFP based PC expansion, ideallya combination of these methods could be used simultaneously to further reduce the computationalcost of PC expansion.
12 Conclusion
We have proposed a procedure for constructing a PC expansion using weighted approximateFekete points and analyzed a model describing oxyhemoglobin dissociation using this procedure.The procedure reduces the required sampling factor for performing PC expansion from multiplicativeto additive. Results from the particular physiological model used in this paper have shown that,when based on weighted approximate Fekete points, a PC expansion using an additive samplingfactor performs comparably to one using a multiplicative factor. Such a result suggests that amoderate amount of data is enough to construct a meaningful PC expansion for prediction and UQanalysis. This procedure is also useful for comparing existing models to clinical observation andexpert behavior in order to analyze model fidelity. ppendix A: Benchmark Testing There are several test functions which are useful for validating a particular implementation of PCexpansion because their sensitivity values can be derived analytically.[57] The functions consideredhere are a polynomial model, the Ishigami function, and the G-function. Since these functions haveanalytical sensitivity values, they are useful for benchmark testing implementation of PC expansioncomputer methods. Sudret et al. have reported the analytical sensitivity values for these functionspreviously.[57] Here we benchmark test a PC expansion as proposed in this paper to confirm the PCexpansion performs adequately. Relative error for estimates of the sensitivity indices was calculatedusing (16).
A.1 Polynomial model
First, we consider the model: Y = 12 n n (cid:89) i =1 (3 X i + 1) , (A.1)for n = 3 where the input variables X i are uniformly distributed over [0, 1]. A polynomial model isuseful for benchmark testing to determine whether the PC expansion procedure can produce exactresults. Given the model in (A.1) is a polynomial of degree 6, a maximal order of p = 6 is chosen.When using the procedure outlined in this paper, including an additive sampling factor f AS = 10,the exact values for the sensitivities were found (tested to 12-digit accuracy).One challenge with PC methods is the choice of polynomial order p . Non-adaptive and adaptivemethods are common; in this manuscript we choose non-adaptive methods for simplicity, where p ischosen a priori . Adaptive methods that refine the polynomial approximation space are more flexiblesince in principle they can be used to drive PC error down to a desired tolerance; however, bothcomputational strategies and convergence for adaptive methods is still an active area of research.To gain an understanding of how the choice of polynomial order p can affect simulation accuracy,we test expansions with p = 3 , , p = 3 isable to estimate the main sensitivities with a relative error less than 1% while the second-order andtotal sensitivity relative errors were within 5% (absolute error within ± . p = 4 estimated the main and second-order sensitivities with a relative error within1% and the largest total sensitivity error of 9.00% (absolute error within ± . p = 5, the PC expansion estimated all sensitivities with relative error within 0.5%. A.2 Ishigami function
Next, we consider the Ishigami function: Y = sin X + a sin X + bX sin X , (A.2)where input variables X i are uniformly distributed over [- π , π ] with a = 7 and b = 0 .
1. All threeinput variables of the Ishigami function are influential (whether directly or through interactions).Although this is uncommon in biomedical engineering applications, this characteristic is useful forbenchmark testing because it typically requires a high order of maximal polynomial degree. TheIshigami function is also commonly used for benchmark testing because of its strong nonlinearityand nonmonoticity.For this function, we tested expansions with maximal order p = { , , , } (see Table 6). Witha threshold of 0.06 for identifying influential sensitivities, a maximal order p = 5 is required toidentify influential parameters and parameter interactions correctly (maximum relative error of40%). A maximal order p = 7 obtains all sensitivities with a relative error within 10%. Increasingthe maximal order to p = 9 estimates all sensitivities with a relative error within 1%. .3 G-function Lastly, we consider Sobol’s G-function: Y = q (cid:89) i =1 | X i − | + a i a i , (A.3)where input variables X i are uniformly distributed over [0, 1], q = 8, and a = { , , , , , , , } .With q = 8, this model can be used to first screen for influential parameters using a low maximalorder PC expansion after which the dimension of the model can be reduced by fixing noninfluentialparameters to their mean value (0.5 in this case). Reducing the dimension lowers the computationalcost for constructing a higher maximal order PC expansion focused on estimating sensitivities forinfluential parameters only.For the G-function, we screened for influential parameters using a maximal order p = 2 whichrequired 55 model runs. Threshold values of 0.02 for main sensitivity and 0.07 for total sensitivitywere required to correctly distinguish between influential ( X − X ) and noninfluential ( X − X )parameters (see Table 7). As planned, This benchmark test then proceeded by setting X − X to their mean value of 0.5, effectively reducing the dimensions of the model, after which a PCexpansion representing the reduced model was constructed to further investigate the sensitivity ofthe model output to input parameters X − X .For the reduced G-function, we tested expansions with maximal order p = { , , , } (Table8). A maximal order of p = 5 is required to rank influential sensitivities correctly. Increasing themaximal order to p = 9 estimated all main and second-order sensitivities with an absolute errorwithin 0.02 and total sensitivities with an absolute error within 0.04. A.4 Conclusion
We have used benchmark testing to analyze the quality of the procedure for constructing a PCexpansion as described in this paper. This PC expansion, which uses an additive sampling factor,performed comparably to other PC expansions previously tested and reported in the literature.Therefore, we conclude that it is prepared for performing analysis of additional models includingthe model for describing oxyhemoglobin dissociation.
13 AcknowledgementsReferences [1]
T. S. I. at the University of Utah , UncertainSCI , July 2020.[2]
C. Bauer and E. Schr¨oder , Carbamino compounds of haemoglobin in human adult andfoetal blood , J Physiol, 227 (1972), pp. 457–471.[3]
R. Benesch and R. Benesch , The effect of organic phosphates from the human erythrocyteon the oxygen equilibrium of human erythrocytes , Biochem Biophys Res Commun, 26 (1967),pp. 162–167.[4]
G. Blatman , Adaptive sparse polynomial chaos expansions for uncertainty propagation andsensitivity analysis , PhD thesis, Clermont-Ferrand 2, 2009.[5]
G. Blatman and B. Sudret , An adaptive algorithm to build up sparse polynomial chaosexpansions for stochastic finite element analysis , Probabilist Eng Mech, 25 (2010), pp. 183–197.[6] ,
Efficient computation of global sensitivity indices using sparse polynomial chaos expan-sions , Reliab Eng Syst Safe, 95 (2010), pp. 1216–1229.[7]
T. Bloom, L. Bos, and N. Levenberg , The Asymptotics of Optimal Designs for PolynomialRegression , arXiv:1112.3735 [stat], (2011). arXiv: 1112.3735.[8]
T. Bloom, L. Bos, N. Levenberg, and S. Waldron , On the Convergence of OptimalMeasures , Constructive Approximation, 32 (2010), pp. 159–179. L. Bos, S. De Marchi, A. Sommariva, and M. Vianello , Computing multivariate Feketeand Leja Points by numerical linear algebra , SIAM J Numer Anal, 48 (2010), p. 1984.[10]
L. M. M. v. d. Bos, B. Sanderse, W. A. A. M. Bierbooms, and G. J. W. van Bus-sel , Generating nested quadrature rules with positive weights based on arbitrary sample sets ,arXiv:1809.09842 [math], (2018). arXiv: 1809.09842.[11]
D. G. Buerk and E. W. Bridges , A simplified algorithm for computing the variation in oxy-hemoglobin saturation with pH, PCO2, T and DPG , Chem Eng Commun, 47 (1986), pp. 113–124.[12]
E. J. Cands, J. K. Romberg, and T. Tao , Stable signal recovery from incomplete and inaccu-rate measurements , Communications on Pure and Applied Mathematics, 59 (2006), pp. 1207–1223.[13]
A. Chanutin and R. R. Curnish , Effect of organic and inorganic phosphates on the oxygenequilibrium of human erythrocytes , Arch Biochem Biophys, 121 (1967), pp. 96–102.[14]
A. Cohen, M. A. Davenport, and D. Leviatan , On the Stability and Accuracy of LeastSquares Approximations , Foundations of Computational Mathematics, 13 (2013), pp. 819–834.[15]
A. Cohen and G. Migliorati , Optimal weighted least-squares methods , SMAI J ComputMath, 3 (2017), pp. 181–203. arxiv:1608.00512 [math.NA].[16]
T. Crestaux, O. Le Maitre, and J.-M. Martinez , Polynomial chaos expansion for sensi-tivity analysis , Reliab Eng Syst Safe, 94 (2009), pp. 1161–1172.[17]
R. K. Dash and J. B. Bassingthwaighte , Erratum to: blood HbO2 and HbCO2 dissociationcurves at varied O2, CO2, pH, 2,3-DPG and temperature levels , Ann Biomed Eng, 38 (2010),pp. 1683–1701.[18]
R. K. Dash, B. Korman, and J. B. Bassingthwaighte , Simple accurate mathematical mod-els of blood HbO2 and HbCO2 dissociation curves at varied physiological conditions: Evaluationand comparison with other models , Eur J Appl Physiol, 116 (2016), pp. 97–113.[19]
W. Donders, W. Huberts, F. van de Vosse, and T. Delhaas , Personalization of modelswith many model parameters: an efficient sensitivity analysis approach , Int J Numer MethodBiomed Eng, 31 (2015).[20]
N. R. Draper and H. Smith , Applied Regression Analysis, 3rd Edition , Wiley, third edi-tion ed., 1998.[21]
S. Dubreuil, M. Berveiller, F. Petitjean, and M. Sala¨un , Construction of bootstrapconfidence intervals on sensitivity indices computed by polynomial chaos expansion , Reliab EngSyst Safe, 121 (2014), pp. 263–275.[22]
V. G. Eck, W. P. Donders, J. Sturdy, J. Feinberg, T. Delhaas, L. R. Hellevik, andW. Huberts , A guide to uncertainty quantification and sensitivity analysis for cardiovascularapplications , Int J Numer Method Biomed Eng, 32 (2016), p. e02755.[23]
S. E. Geneser, R. M. Kirby, and R. S. MacLeod , Application of stochastic finite elementmethods to study the sensitivity of ECG forward modeling to organ conductivity , IEEE TransBiomed Eng, 55 (2008), pp. 31–40.[24]
L. Guo, A. Narayan, L. Yan, and T. Zhou , Weighted approximate Fekete points: Samplingfor least-squares polynomial approximation , SIAM J Sci Comput, 40 (2018), pp. A366–A387.arXiv:1708.01296 [math.NA].[25]
M. Hadigol and A. Doostan , Least squares polynomial chaos expansion: A review ofsampling strategies , Computer Methods in Applied Mechanics and Engineering, 332 (2018),pp. 382–407.[26]
H. C. Hemmings and P. M. Hopkins , Foundations of anesthesia: basic sciences for clinicalpractice , Elsevier Health Sciences, 2006.[27]
M. P. Hlastala, R. D. Woodson, and B. Wranne , Influence of temperature on hemoglobin-ligand interaction in whole blood , J Appl Physiol, 43 (1977), pp. 545–550. B. Iooss, A. Janon, G. Pujol, with contributions from Khalid Boumhaout, S. D.Veiga, T. Delage, J. Fruth, L. Gilquin, J. Guillaume, L. Le Gratiet, P. Lemaitre,B. L. Nelson, F. Monari, R. Oomen, O. Rakovec, B. Ramos, O. Roustant, E. Song,J. Staum, R. Sueur, T. Touati, and F. Weber , Sensitivity: Global Sensitivity Analysis ofModel Outputs , 2018. R package version 1.15.2.[29]
N. Joels and L. Pugh , The carbon monoxide dissociation curve of human blood , J Physiol,142 (1958), pp. 63–77.[30]
L. R. Johnson , Essential medical physiology , Academic Press, 2003.[31]
G. R. Kelman , Digital computer subroutine for the conversion of oxygen tension into satura-tion , J Appl Physiol, 21 (1966), pp. 1375–1376.[32]
M. C. Kennedy and A. O’Hagan , Bayesian calibration of computer models , Journal of theRoyal Statistical Society: Series B (Statistical Methodology), 63 (2001), pp. 425–464.[33]
V. Keshavarzzadeh, R. Kirby, and A. Narayan , Numerical Integration in Multiple Dimen-sions with Designed Quadrature , SIAM Journal on Scientific Computing, 40 (2018), pp. A2033–A2061. arXiv:1804.06501 [cs.NA].[34]
M. V. L Bos, F. Piazzon , Near optimal polynomial regression on norming meshes . SampTA2019, IEEE Xplore Digital Library, to appear.[35]
M. Matej´ak, T. Kulh´anek, and S. Matouˇsek , Adair-based hemoglobin equilibrium withoxygen, carbon dioxide and hydrogen ion activity , Scand J Clin Lab Invest, 75 (2015), pp. 113–120.[36]
J. Matthew, J. Morrow, R. Wittebort, and F. Gurd , Quantitative determination ofcarbamino adducts of alpha and beta chains in human adult hemoglobin in presence and absenceof carbon monoxide and 2,3-diphosphoglycerate , J Biol Chem, 252 (1977), pp. 2234–2244.[37]
S. Mehrotra and D. Papp , Generating Moment Matching Scenarios Using OptimizationTechniques , SIAM Journal on Optimization, 23 (2013), pp. 963–999.[38]
N. Naeraa, E. S. Petersen, and E. Boye , The influence of simultaneous, independentchanges in pH and carbon dioxide tension on the in vitro oxygen tension-saturation relationshipof human blood , Scand J Clin Lab Invest, 15 (1963), pp. 141–151.[39]
A. Narayan , Computation of induced orthogonal polynomial distributions , Electron TransNumer Anal, 50 (2018), pp. 71–97. arXiv:1704.08465 [math].[40]
A. Narayan, J. Jakeman, and T. Zhou , A Christoffel function weighted least squares algo-rithm for collocation approximations , Mathematics of Computation, 86 (2017), pp. 1913–1947.arXiv: 1412.4305 [math.NA].[41]
A. Narayan and D. Xiu , Constructing Nested Nodal Sets for Multivariate Polynomial Inter-polation , SIAM Journal on Scientific Computing, 35 (2013), pp. A2293–A2315.[42]
A. Narayan and T. Zhou , Stochastic Collocation on Unstructured Multivariate Meshes , Com-munications in Computational Physics, 18 (2015), pp. 1–36. arXiv:1501.05891 [math.NA].[43]
J. Oakley and A. O’hagan , Bayesian inference for the uncertainty distribution of computermodel outputs , Biometrika, 89 (2002), pp. 769–784.[44]
H. Oudemans-van Straaten, G. Scheffer, and C. Stoutenbeek , Analysis of P50 andoxygen transport in patients after cardiac surgery , Intensive Care Med, 22 (1996), pp. 781–789.[45]
F. Piazzon, A. Sommariva, and M. Vianello , Caratheodory-tchakaloff least squares , in2017 International Conference on Sampling Theory and Applications (SampTA), July 2017,pp. 672–676.[46]
H. Rauhut and S. Foucart , A Mathematical Introduction to Compressive Sensing ,Birkhuser, New York, 2013 edition ed., Aug. 2013.[47]
R. B. Reeves , The effect of temperature on the oxygen equilibrium curve of human blood ,Respir Physiol, 42 (1980), pp. 317–328.[48]
R. A. Rhoades and D. R. Bell , Medical physiology: Principles for clinical medicine , Lip-pincott Williams & Wilkins, 2012. T. Rosenthal , The effect of temperature on the pH of blood and plasma in vitro , J Biol Chem,173 (1948), pp. 25–30.[50]
A. Saltelli , Making best use of model evaluations to compute sensitivity indices , ComputPhys Commun, 145 (2002), pp. 280–297.[51]
J. Severinghaus , Oxyhemoglobin dissociation curve correction for temperature and pH vari-ation in human blood , J Appl Physiol, 12 (1958), pp. 485–486.[52]
J. W. Severinghaus , Simple, accurate equations for human blood O2 dissociation computa-tions , J Appl Physiol, 46 (1979), pp. 599–602.[53]
Y. Shin and D. Xiu , Nonadaptive Quasi-Optimal Points Selection for Least Squares LinearRegression , SIAM Journal on Scientific Computing, 38 (2016), pp. A385–A411.[54]
O. Siggaard-Andersen , Oxygen-linked hydrogen ion binding of human hemoglobin. Effects ofcarbon dioxide and 2,3-diphosphoglycerate I. studies on erythrolysate , Scand J Clin Lab Invest,27 (1971), pp. 351–360.[55]
O. Siggaard-Andersen, P. Wimberley, I. G¨othgen, and M. Siggaard-Andersen , Amathematical model of the hemoglobin-oxygen dissociation curve of human blood and of theoxygen partial pressure as a function of temperature , Clin Chem, 30 (1984), pp. 1646–1651.[56]
I. Sobol , Global sensitivity indices for nonlinear mathematical models and their Monte Carloestimates , Math Comput Simulat, 55 (2001), pp. 271–280.[57]
B. Sudret , Global sensitivity analysis using polynomial chaos expansions , Reliab Eng SystSafe, 93 (2008), pp. 964–979.[58]
H. M. Thomas, S. S. Lefrak, R. S. Irwin, H. W. Fritts, and P. R. Caldwell , Theoxyhemoglobin dissociation curve in health and disease: Role of 2,3-diphosphoglycerate , Am JMed, 57 (1974), pp. 331–348.[59]
I. Tyuma and K. Shimizu , Different response to organic phosphates of human fetal and adulthemoglobins. , Arch Biochem Biophys, 129 (1969), pp. 404–405.[60]
R. M. Winslow, M. Samaja, N. J. Winslow, L. Rossi-Bernardi, and R. I. Shrager , Simulation of continuous blood O2 equilibrium curve over physiological pH, DPG, and PCO2range , J Appl Physiol, 54 (1983), pp. 524–529.[61]
D. Xiu , Numerical methods for stochastic computations: A spectral method approach , PrincetonUniversity Press, July 2010.[62]
D. Xiu and G. E. Karniadakis , The Wiener–Askey polynomial chaos for stochastic differ-ential equations , SIAM J Sci Comput, 24 (2002), pp. 619–644., SIAM J Sci Comput, 24 (2002), pp. 619–644.