# Efficient Basis Change for Sparse-Grid Interpolating Polynomials with Application to T-Cell Sensitivity Analysis

Authors: Gregery T. Buzzard
Journal: Computational Biology Journal
 (2013)

---

## Abstract

Sparse-grid interpolation provides good approximations to smooth functions in high dimensions based on relatively few function evaluations, but in standard form it is expressed in Lagrange polynomials. Here, we give a block-diagonal factorization of the basis-change matrix to give an efficient conversion of a sparse-grid interpolant to a tensored orthogonal polynomial (or gPC) representation. We describe how to use this representation to give an efficient method for estimating Sobol' sensitivity coefficients and apply this method to analyze and efficiently approximate a complex model of T-cell signaling events.


---

## Body

## 1. Introduction

A common problem in many areas of computational mathematics is to approximate a given function based on a small number of functional evaluations or observations. This problem arises in numerical methods for PDE [1, 2], sensitivity analysis [3–5], uncertainty quantification [2], many areas of modeling [6, 7], and other settings. As a result, there are a large number of approaches to this problem, and the literature is large and growing quickly. In settings in which the points of evaluation may be chosen at will, two common approaches are generalized polynomial chaos (gPC) using cubature and sparse grid collocation [8]. In other settings in which the points of evaluation are given, common approaches include RS-HDMR, cut-HDMR, ANOVA decomposition, kriging, and moving least squares; see, for example, [9, Chapter 5] for a discussion of such methods and further references.Sparse grid collocation has been used widely in recent years as a means of providing a reasonable approximation to a smooth function,f, defined on a hypercube in ℝn, based on relatively few function evaluations [2]. This method produces a polynomial interpolant using Lagrange interpolating polynomials based on function values at points in a union of product grids of small dimension [10, 11]. Using barycentric interpolation to evaluate the resulting polynomial [12], this method is a viable alternative to an expansion of f in terms of a sum of products of one-dimensional orthogonal polynomials. This latter approach is known as generalized polynomial chaos (gPC) or spectral decomposition and is obtained via standard weighted L2 techniques. However, the orthogonality implicit in the gPC representation often provides many advantages over the Lagrange representation, particularly in applications to differential equations, in which the gPC representation is closely related to spectral methods. Other advantages of the gPC representation include the ability to estimate convergence as more points are added to the sparse grid, and the ability to estimate variance-based sensitivity coefficients quickly and accurately [4]. A common approach to obtain the gPC coefficients is to use numerical integration by applying a cubature rule. However, cubature rules to integrate the product of two polynomials up to the degree in the sparse grid interpolant typically require either more or different points than found in the sparse grid itself.In this paper, we provide an efficient algorithm for converting from the Lagrange interpolating polynomial to an equivalent gPC polynomial using only the function values at the sparse grid points. The foundation of this algorithm is a matrix factorization based on the fact that the sparse grid is a union of small product grids. More precisely, letΦ be the matrix obtained by evaluating each of the gPC basis functions (one per column) at each of the sparse grid points (one per row). This matrix produces the gPC coefficients, c, from the function values, f, by solving Φc=f for c. We show below that Φ-1 factors into a product of block diagonal matrices in which each block corresponds to one of the small product grids composing the full sparse grid. We also demonstrate the utility of this method for calculating variance-based sensitivity coefficients on a model of T-cell signaling.Methods for changing basis between sets of orthogonal polynomials are described in many places, including [13–15]. Most of these results focus on the case of polynomials of one variable, and none focuses on the case of polynomial bases for sparse grid interpolants. For further background, Boyd [16, Chapter 5] discusses the matrix multiplication transform for converting between Lagrange and spectral representations, while Chapter 10 discusses the special case of the Chebyshev polynomial expansion using the FFT.Organization. In Section 2, we provide some background into Smolyak’s algorithm for sparse grid interpolation and the calculation of Sobol’ sensitivity coefficients. In Section 3.1, we give the factorization of Φ-1 and in Section 3.2 we describe the algorithm to produce the gPC representation based on function values at the sparse grid points and give an upper bound on the complexity of the algorithm. In Section 4.1, we give numerical results on running time and accuracy of the algorithm for a simple test function. In Section 4.2, we describe the application of this algorithm to compute Sobol’ sensitivity coefficients for a model of T-cell signaling.


## 2. Background

### 2.1. Sparse Grid Interpolation

In this section, we provide some background on Smolyak’s algorithm for sparse grid interpolation. The discussion here is based on [10]. The foundation for sparse grid interpolation is interpolation in one dimension using Lagrange interpolation as follows:
(1)Ui(f)=∑j=1mi‍f(xji)Lji,
where i∈ℕ,  xji∈[-1,1], and Lji is the Lagrange polynomial satisfying Lji(xki)=δjk. A common choice for sparse grid interpolation is to use the Chebyshev-Gauss-Lobatto (CGL) points, in which case mi=2i-1+1 and xji=-cos((j-1)π/2i-1) for i>1, since this choice provides a nesting property that is vital for efficient interpolation in higher dimensions. Moreover, m1=1 and x11=0 for this choice. Other choices, including Gauss-Patterson nodes and nodes on a (semi)infinite interval as in [17], are also possible; the methods below apply for these choices as well. These one-dimensional formulas may be tensored in dimension d>1 to yield
(2)Ui(f)=(Ui1⊗⋯⊗Uid)(f)=∑1≤j≤mi‍f(xji)Lji,
where i and j are multi-indices with componentwise partial order, 1 is the multi-index of all 1s, xji=(xj1i1,…,xjdid), and Lji(x)=Lj1i1(x1)⋯Ljdid(xd). This formula requires mi1⋯mid function values sampled on a product grid. Note, however, that when some of the ik are 1, then this grid is of dimension less than d (since m1=1).Linear combinations of these formulas produce the Smolyak formulas. LetU0=0 and Δi=Ui-Ui-1 for i∈ℕ, and define |i|=i1+⋯+id. Then for q≥d, we have
(3)A(q,d)=∑d≤|i|≤q‍Δi,
where Δi is the tensor product of the Δik. A multinomial expansion and q=d+k produce
(4)A(d+k,d)=∑k+1≤|i|≤d+k‍(-1)d+k-|i|(d-1|i|-k-1)Ui.
This formula is known as the combination technique, which first appeared in connection with sparse grids in [18]. When d>k, as is common for large d, we can replace k+1≤|i| by d≤|i|. An anisotropic version of this formula is described in [19]. The formula for this version is largely the same as (3) with the index set d≤|i|≤q replaced by a more general index set I, characterized by the property that if i∈I and ik>1, then i-ek∈I, where ek has 1 in the kth position and zeros elsewhere.Important results from [10] include that A(d+k,d)(P)=P for all polynomials of degree at most k, and that if the derivatives Dαf are continuous for all multi-indices α with αi≤m for all i, then
(5)∥f-A(q,d)(f)∥∞≤cd,mnm(logn)(m+2)(d-1)+1∥f∥∞,
where n=n(q,d) is the number of points in the sparse grid, X(q,d), for A(q,d). There is a similar error estimate for a weighted L2 norm of the Chebyshev expansion of f-A(q,d)(f). They note also that n(d+k,d) is asymptotic to 2kdk/k! as d tends to ∞, in the sense that the ratio tends to 1. More precise estimates on convergence may be found in [20].


### 2.2. Global Sensitivity Analysis

A common problem in analyzing complex biological models with unknown parameters is identifying which parameters play an important role in modulating the response of the model. The variance-based sensitivity coefficients of Sobol’ [21] can provide a great deal of insight into this problem.The calculation of these coefficients for a functionf(x)=f(x1,…,xn) defined on a hypercube in ℝn starts by writing f as a normalized sum of functions that depend on a specified subset of the following variables:
(6)f(x)=f0+∑i=1n‍fi(xi)+∑1≤i<j≤n‍fij(xi,xj)+⋯+f12⋯n(x1,…,xn),
with a weighted orthogonality condition imposed on pairs of functions in this decomposition and a 0 mean condition in each variable separately imposed on each function individually. Depending on the context, this decomposition is known as the Sobol’ or ANOVA decomposition of f or the HDMR representation of f. One method for obtaining this decomposition is by calculating the gPC expansion of f [4, 7]. In fact, the gPC expansion is essentially a refinement of the above decomposition in which each component function in the sum above is written as a sum of a product of one-dimensional orthogonal polynomials.The Sobol’ or variance-based sensitivity coefficients forf are then defined in terms of the variance of the functions in the Sobol’ decomposition. More precisely, let Hd denote the hypercube [-1,1]d and assume that the decomposition above takes place on this set. Let Vd=2d be the volume of Hd, and for a function g on this set, define the variance of g to be
(7)D(g)=1Vd∫Hd‍g2(x)dx-(1Vd∫Hd‍g(x)dx)2.
In terms of the Sobol’ decomposition given above, the main effect sensitivity coefficient for xj is Sj=D(fj)/D(f). This is generalized to interaction coefficients, SM by replacing fj by fM, where M is any subset of {1,…,n}. Also, the total effect sensitivity coefficient STj is the sum of SM over all M containing j.As shown in [4], the gPC expansion of f leads to an efficient method for calculating these sensitivity coefficients: D(fM) is simply the sum of squares of coefficients of the gPC polynomials in the expansion of fM. Hence, the change of basis algorithm described above provides an efficient way to estimate Sobol’ sensitivity coefficients based on function values at sparse grid nodes. We simply evaluate the function at sparse grid nodes, change to a gPC basis using Legendre polynomials, and find the appropriate sum of squares of coefficients to calculate the Sj and STj.


## 2.1. Sparse Grid Interpolation

In this section, we provide some background on Smolyak’s algorithm for sparse grid interpolation. The discussion here is based on [10]. The foundation for sparse grid interpolation is interpolation in one dimension using Lagrange interpolation as follows:
(1)Ui(f)=∑j=1mi‍f(xji)Lji,
where i∈ℕ,  xji∈[-1,1], and Lji is the Lagrange polynomial satisfying Lji(xki)=δjk. A common choice for sparse grid interpolation is to use the Chebyshev-Gauss-Lobatto (CGL) points, in which case mi=2i-1+1 and xji=-cos((j-1)π/2i-1) for i>1, since this choice provides a nesting property that is vital for efficient interpolation in higher dimensions. Moreover, m1=1 and x11=0 for this choice. Other choices, including Gauss-Patterson nodes and nodes on a (semi)infinite interval as in [17], are also possible; the methods below apply for these choices as well. These one-dimensional formulas may be tensored in dimension d>1 to yield
(2)Ui(f)=(Ui1⊗⋯⊗Uid)(f)=∑1≤j≤mi‍f(xji)Lji,
where i and j are multi-indices with componentwise partial order, 1 is the multi-index of all 1s, xji=(xj1i1,…,xjdid), and Lji(x)=Lj1i1(x1)⋯Ljdid(xd). This formula requires mi1⋯mid function values sampled on a product grid. Note, however, that when some of the ik are 1, then this grid is of dimension less than d (since m1=1).Linear combinations of these formulas produce the Smolyak formulas. LetU0=0 and Δi=Ui-Ui-1 for i∈ℕ, and define |i|=i1+⋯+id. Then for q≥d, we have
(3)A(q,d)=∑d≤|i|≤q‍Δi,
where Δi is the tensor product of the Δik. A multinomial expansion and q=d+k produce
(4)A(d+k,d)=∑k+1≤|i|≤d+k‍(-1)d+k-|i|(d-1|i|-k-1)Ui.
This formula is known as the combination technique, which first appeared in connection with sparse grids in [18]. When d>k, as is common for large d, we can replace k+1≤|i| by d≤|i|. An anisotropic version of this formula is described in [19]. The formula for this version is largely the same as (3) with the index set d≤|i|≤q replaced by a more general index set I, characterized by the property that if i∈I and ik>1, then i-ek∈I, where ek has 1 in the kth position and zeros elsewhere.Important results from [10] include that A(d+k,d)(P)=P for all polynomials of degree at most k, and that if the derivatives Dαf are continuous for all multi-indices α with αi≤m for all i, then
(5)∥f-A(q,d)(f)∥∞≤cd,mnm(logn)(m+2)(d-1)+1∥f∥∞,
where n=n(q,d) is the number of points in the sparse grid, X(q,d), for A(q,d). There is a similar error estimate for a weighted L2 norm of the Chebyshev expansion of f-A(q,d)(f). They note also that n(d+k,d) is asymptotic to 2kdk/k! as d tends to ∞, in the sense that the ratio tends to 1. More precise estimates on convergence may be found in [20].


## 2.2. Global Sensitivity Analysis

A common problem in analyzing complex biological models with unknown parameters is identifying which parameters play an important role in modulating the response of the model. The variance-based sensitivity coefficients of Sobol’ [21] can provide a great deal of insight into this problem.The calculation of these coefficients for a functionf(x)=f(x1,…,xn) defined on a hypercube in ℝn starts by writing f as a normalized sum of functions that depend on a specified subset of the following variables:
(6)f(x)=f0+∑i=1n‍fi(xi)+∑1≤i<j≤n‍fij(xi,xj)+⋯+f12⋯n(x1,…,xn),
with a weighted orthogonality condition imposed on pairs of functions in this decomposition and a 0 mean condition in each variable separately imposed on each function individually. Depending on the context, this decomposition is known as the Sobol’ or ANOVA decomposition of f or the HDMR representation of f. One method for obtaining this decomposition is by calculating the gPC expansion of f [4, 7]. In fact, the gPC expansion is essentially a refinement of the above decomposition in which each component function in the sum above is written as a sum of a product of one-dimensional orthogonal polynomials.The Sobol’ or variance-based sensitivity coefficients forf are then defined in terms of the variance of the functions in the Sobol’ decomposition. More precisely, let Hd denote the hypercube [-1,1]d and assume that the decomposition above takes place on this set. Let Vd=2d be the volume of Hd, and for a function g on this set, define the variance of g to be
(7)D(g)=1Vd∫Hd‍g2(x)dx-(1Vd∫Hd‍g(x)dx)2.
In terms of the Sobol’ decomposition given above, the main effect sensitivity coefficient for xj is Sj=D(fj)/D(f). This is generalized to interaction coefficients, SM by replacing fj by fM, where M is any subset of {1,…,n}. Also, the total effect sensitivity coefficient STj is the sum of SM over all M containing j.As shown in [4], the gPC expansion of f leads to an efficient method for calculating these sensitivity coefficients: D(fM) is simply the sum of squares of coefficients of the gPC polynomials in the expansion of fM. Hence, the change of basis algorithm described above provides an efficient way to estimate Sobol’ sensitivity coefficients based on function values at sparse grid nodes. We simply evaluate the function at sparse grid nodes, change to a gPC basis using Legendre polynomials, and find the appropriate sum of squares of coefficients to calculate the Sj and STj.


## 3. Main Results

### 3.1. Change of Basis Matrix and Factorization

In this section, we describe the matrix factorization that is central to the efficient change to a gPC basis.LetPj, j≥0 be a set of polynomials, each of degree j, orthogonal on the interval [-1,1] using inner product ∫-11‍fgwdx, where w is a nonnegative weight function that is positive in the interior of the interval (straightforward modifications allow for unbounded intervals, but for clarity we restrict to the finite interval).The expansion in (4) together with (1) gives a decomposition of A(q,d) into a sum of products of Lagrange polynomials, with each Ui corresponding to a full product grid, Xi (although usually Xi is less than d-dimensional since each entry of i, that is 1, produces only a single point for the corresponding dimension). On this grid, each nontrivial coordinate direction k has i=ik>1, and the number of points in this direction is mi=2i-1+1. Hence, the degree of each Lji is mi-1, so we write each Lji as a linear combination of P0,…,Pmi-1. For a fixed value of i, we have
(8)Lji=∑m=0mi-1Cj,mi‍Pm.
Taking tensor products as in (1), we obtain a basis in terms of tensor products of Pn for interpolation on the product grid Xi. Taking a union over the Xi that appear implicitly in (4) via Ui, we obtain a basis for interpolation on the entire sparse grid. The problem then is to determine the (gPC) coefficients in this basis for specified function values, f, at the sparse grid nodes.Conceptually, the simplest approach to finding the gPC coefficients starts by evaluating each basis function,P, at each point in X=X(q,d). These values form a matrix, Φ, with each basis function corresponding to one column and each point in X corresponding to one row.Definition 1.
LetPj,  j=1,…,n be the gPC basis polynomials for the sparse grid, X, having points x1,…,xn. Then the matrix Φ is the matrix with entries Φi,j=Pj(xi).With this definition, and a gPC coefficient vector,c, the product Φc gives the value at each point in X of the polynomial determined by the coefficients in c. Hence, one can solve the interpolation problem by solving Φc=f for c. Of course, simply constructing this matrix takes O(n2) operations, while solving the system in this form using elimination takes O(n3) operations.The alternative presented here is to give a factorization ofΦ-1 based on (4).Theorem 2.
The change-of-basis matrixΦ-1 has a factorization
(9)Φ-1=πTD^wΦ^-1π,
where Φ^ is block-diagonal, with each block corresponding to a subgrid, Xi; D^w is diagonal; π is a matrix that includes ℝn into ℝN for some N>n.Proof.
To start, we order the columns ofΦ so that for each multi-index, i, the set of column indices of basis functions for interpolation on Xi from (8) are the same as the set of row indices for points in Xi. Let n be the number of points in X and ni the number in Xi. Let πi be the projection from ℝn to ℝni obtained by restricting to the indices corresponding to Xi. In this case, Φi=πiΦπiT is the matrix obtained by evaluating basis functions at points, with the basis functions and points restricted to those associated with Xi. Hence, given f=f(X), we can interpolate f on Xi by solving Φici=πif for ci.
SinceUi(f) in (1) also interpolates f on Xi, we have
(10)Ui(f)(Xi)=Φici=πiΦπiTci.
Since the gPC basis functions associated with Xi were obtained by changing basis from the Lagrange representation of Ui, we see that Ui(f) and the polynomial with gPC coefficients πiTci are the same polynomial. Hence, we may extend the equality in (10) to all of X by replacing Xi by X on the left and dropping πi on the right to obtain
(11)Ui(f)(X)=ΦπiTci.
Let wi denote the weight for Ui in (4). Using (11) with (4), we have
(12)f=A(f)(X)=∑i∈I‍wiUi(f)(X)=∑i‍wiΦπiTci=Φ(∑i‍wiπiTci).
Recalling that ci=Φi-1πif, we have
(13)Φ-1f=∑i∈I‍wiπiTΦi-1πif.
Let i1,…,im be an enumeration of the indices in I, let π be the matrix obtained by vertical concatenation of πi1,…,πim, let Φ^ be the block diagonal matrix with Φi1,…,Φim as blocks, and let D^w be the diagonal matrix with diagonal entries wi1,…,wim, each repeated according to the size of the corresponding block. Then (13) gives the factorization
(14)Φ-1=πTD^wΦ^-1π.


### 3.2. Algorithm and Complexity

Based on the previous section, the algorithm to find the coefficient vector,c, to represent the sparse-grid interpolating polynomial in gPC basis is simply to use (14) to find c=Φ-1f. We avoid the matrix inverse by defining f^=πf, solving Φ^c^=f^ for c^, and then using c=πTD^wc^.In this section, we show that for fixedk in (4) and large dimension, d, the complexity of this algorithm is linear in the number of points of evaluation. We frame this result as a corollary to the following proposition, which gives a bound on the total number of operations needed to compute on each subgrid, Xi, in turn, as needed to solve Φ^c^=f^ for c^.In both of the following results,n(d+k,d) is the number of points in the sparse grid associated with A(d+k,d).Proposition 3.
Letm≥1 and k≥1 be fixed. Then there is ck,m>0 so that if d>k, then
(15)∑d≤|i|≤d+k‍|Xi|m≤ck,mn(d+k,d).In this proposition, we use the fact mentioned above that whend>k, then the sparse grid X is a union of Xi over i as indicated in the sum given here. Given this result, the following theorem is nearly immediate.Theorem 4.
For fixedk, there is ck>0 so that for d>k, the coefficients in the gPC expansion of A(d+k,d)(f) may be found using (14) with computation time bounded by ckn(d+k,d); that is, for fixed k, the running time is linear in the number of grid points.We have not attempted to calculate the best possibleck in the proof of this result. In fact, numerical results given below show that in some cases, the coefficient ck actually decreases as k increases due to the spread of fixed overhead costs over a small number of points when k and d are small. Additionally, the numerical results show that the time per point evaluated is roughly constant through k=8.Note that this algorithm computes the gPC coefficients of the interpolating polynomial rather than the original function itself. However, for aCr function, f, with r≥1, the error estimate from (5) implies that the interpolating function converges uniformly to f as the depth, k, increases. Since the gPC coefficients are obtained by weighted integration of f against the orthogonal polynomials, the gPC coefficients for the interpolating polynomials converge to the gPC coefficients for f, and (5) provides a means to estimate the error in these coefficients based on the set of orthogonal polynomials and their corresponding weights.Proof of Proposition3.
Given|i|=d+κ between d and d+k, consider the number of points in the corresponding grid Xi. Let r=r(i) be the number of nontrivial entries in i (i.e., entries larger than 1). Suppose without loss that the first r entries of i are nontrivial with entries ij=κj+1≥2. Then, the number of points in this grid is ∏j=1r‍(2κj+1). If κ1 is the only entry larger than 1, then κ1=κ-(r-1), and so the number of points is exactly 3r-1(2κ-r+1+1). A simple counting argument implies that any other distribution of nontrivial entries produces no more grid points, so
(16)|Xi|≤3r-12κ-r+2.
To constructi with |i|=d+κ and r nontrivial entries, we select r out of d nontrivial entries, then distribute κ elements to these entries and so that each entry gets at least one element. There are (dr)(κ-1r-1) ways to do this. Summing from r=1 to κ, and using (dr)≤dr, we have
(17)∑|i|=d+κ‍|Xi|m≤∑r=1κ‍(dr)(κ-1r-1)(3m)r-1(2m)κ-r+2≤(2m)κ+1d∑q=0κ-1‍(κ-1q)(3m)qdq(2m)q.
The summation on the right hand side is a binomial expansion of a power of κ-1, and distributing a factor of (2m)κ-1 into this power converts the right hand side to 4md(3md+2m)κ-1. Using this in place of the sum over i and summing over κ gives
(18)4md∑κ=1k‍(3md+2m)κ-1=4md(3md+2m)k-13md+2m-1≤2dk4m(3m+2md)k-1.
As noted above, n(d+k,d) is asymptotically equal to 2kdk/k! for large d in the sense that for fixed k, limd→∞n(d+k,d)/(2kdk/k!)=1. In particular, there is ak>0 so that dk≤akn(d+k,d) for all d≥1. Hence, the results just given imply
(19)∑κ=1k‍∑|i|=d+κ‍|Xi|m≤2akn(d+k,k)4m(2·3m)k-1.
Hence, taking ck,m=2ak4m(2·3m)k-1, we obtain the desired result.Proof of Theorem4.
Given the sparse grid,X, associated with A(d+k,k) and the vector of values, f, we need first to form the matrices π, D^w, and Φ^ as in (13). From the construction of π, it has exactly one nonzero entry per row and is n×N, where n=n(d,k) is asymptotic to 2kdk/k! as d increases [10], and N=∑d≤|i|≤d+k‍|Xi|. Hence, π and the diagonal matrix, D^w, may be constructed and applied in time O(N).
A given block ofΦ^ is obtained by evaluating products of one-dimensional polynomials on the points of a product grid, Xi. For a nontrivial entry ij=1+κj, the polynomials in xj have degree at most 2κj. Using a three-term recurrence, these polynomials may be evaluated at a given point in O(2κj). Since the projection of Xi to the xj coordinate has 2κj+1 points, these polynomials are evaluated in time O(22κj). This is repeated for each nontrivial κj for a total of O(k22κj), which is O(k|Xi|2). Then for each point in Xi, we multiply at most k polynomials, which are O(k|Xi|). Summing over i and applying Proposition 3, we see that Φ^ may be constructed in time O(kn(d+k,d)).
As noted above, to find the gPC coefficients,c, we define f^=πf, solve Φ^c^=f^ for c^, and then use c=πTD^wc^. Since Φ^ has blocks of size |Xi|×|Xi|, we may solve each block in time at most O(|Xi|3). Combining this with Proposition 3 and the estimates O(N) and O(kN) above, we obtain the theorem.A practical point is that many of the blocks inΦ^ are identical up to a permutation of the rows. Hence, for each of these blocks, we may use a single LU decomposition and appropriate permutations of the entries in f to reduce the total running time (although not, of course, the linear dependence described in the theorem).We note also that the algorithm given here is compatible with the anisotropic adaptive sparse grids of [19]. The analysis given here shows that the time per point depends only on the maximum size of a block in the factorization of Φ-1, and this size is bounded by the maximum depth, k, in both the isotropic and the anisotropic cases.Finally, the analysis for fixedk and increasing d is relevant in that as models become more complex, d increases, while often the important behavior of a model can be captured by a relatively low-order approximation. An example of this is seen in the T-cell signaling model in Section 4.2.


## 3.1. Change of Basis Matrix and Factorization

In this section, we describe the matrix factorization that is central to the efficient change to a gPC basis.LetPj, j≥0 be a set of polynomials, each of degree j, orthogonal on the interval [-1,1] using inner product ∫-11‍fgwdx, where w is a nonnegative weight function that is positive in the interior of the interval (straightforward modifications allow for unbounded intervals, but for clarity we restrict to the finite interval).The expansion in (4) together with (1) gives a decomposition of A(q,d) into a sum of products of Lagrange polynomials, with each Ui corresponding to a full product grid, Xi (although usually Xi is less than d-dimensional since each entry of i, that is 1, produces only a single point for the corresponding dimension). On this grid, each nontrivial coordinate direction k has i=ik>1, and the number of points in this direction is mi=2i-1+1. Hence, the degree of each Lji is mi-1, so we write each Lji as a linear combination of P0,…,Pmi-1. For a fixed value of i, we have
(8)Lji=∑m=0mi-1Cj,mi‍Pm.
Taking tensor products as in (1), we obtain a basis in terms of tensor products of Pn for interpolation on the product grid Xi. Taking a union over the Xi that appear implicitly in (4) via Ui, we obtain a basis for interpolation on the entire sparse grid. The problem then is to determine the (gPC) coefficients in this basis for specified function values, f, at the sparse grid nodes.Conceptually, the simplest approach to finding the gPC coefficients starts by evaluating each basis function,P, at each point in X=X(q,d). These values form a matrix, Φ, with each basis function corresponding to one column and each point in X corresponding to one row.Definition 1.
LetPj,  j=1,…,n be the gPC basis polynomials for the sparse grid, X, having points x1,…,xn. Then the matrix Φ is the matrix with entries Φi,j=Pj(xi).With this definition, and a gPC coefficient vector,c, the product Φc gives the value at each point in X of the polynomial determined by the coefficients in c. Hence, one can solve the interpolation problem by solving Φc=f for c. Of course, simply constructing this matrix takes O(n2) operations, while solving the system in this form using elimination takes O(n3) operations.The alternative presented here is to give a factorization ofΦ-1 based on (4).Theorem 2.
The change-of-basis matrixΦ-1 has a factorization
(9)Φ-1=πTD^wΦ^-1π,
where Φ^ is block-diagonal, with each block corresponding to a subgrid, Xi; D^w is diagonal; π is a matrix that includes ℝn into ℝN for some N>n.Proof.
To start, we order the columns ofΦ so that for each multi-index, i, the set of column indices of basis functions for interpolation on Xi from (8) are the same as the set of row indices for points in Xi. Let n be the number of points in X and ni the number in Xi. Let πi be the projection from ℝn to ℝni obtained by restricting to the indices corresponding to Xi. In this case, Φi=πiΦπiT is the matrix obtained by evaluating basis functions at points, with the basis functions and points restricted to those associated with Xi. Hence, given f=f(X), we can interpolate f on Xi by solving Φici=πif for ci.
SinceUi(f) in (1) also interpolates f on Xi, we have
(10)Ui(f)(Xi)=Φici=πiΦπiTci.
Since the gPC basis functions associated with Xi were obtained by changing basis from the Lagrange representation of Ui, we see that Ui(f) and the polynomial with gPC coefficients πiTci are the same polynomial. Hence, we may extend the equality in (10) to all of X by replacing Xi by X on the left and dropping πi on the right to obtain
(11)Ui(f)(X)=ΦπiTci.
Let wi denote the weight for Ui in (4). Using (11) with (4), we have
(12)f=A(f)(X)=∑i∈I‍wiUi(f)(X)=∑i‍wiΦπiTci=Φ(∑i‍wiπiTci).
Recalling that ci=Φi-1πif, we have
(13)Φ-1f=∑i∈I‍wiπiTΦi-1πif.
Let i1,…,im be an enumeration of the indices in I, let π be the matrix obtained by vertical concatenation of πi1,…,πim, let Φ^ be the block diagonal matrix with Φi1,…,Φim as blocks, and let D^w be the diagonal matrix with diagonal entries wi1,…,wim, each repeated according to the size of the corresponding block. Then (13) gives the factorization
(14)Φ-1=πTD^wΦ^-1π.


## 3.2. Algorithm and Complexity

Based on the previous section, the algorithm to find the coefficient vector,c, to represent the sparse-grid interpolating polynomial in gPC basis is simply to use (14) to find c=Φ-1f. We avoid the matrix inverse by defining f^=πf, solving Φ^c^=f^ for c^, and then using c=πTD^wc^.In this section, we show that for fixedk in (4) and large dimension, d, the complexity of this algorithm is linear in the number of points of evaluation. We frame this result as a corollary to the following proposition, which gives a bound on the total number of operations needed to compute on each subgrid, Xi, in turn, as needed to solve Φ^c^=f^ for c^.In both of the following results,n(d+k,d) is the number of points in the sparse grid associated with A(d+k,d).Proposition 3.
Letm≥1 and k≥1 be fixed. Then there is ck,m>0 so that if d>k, then
(15)∑d≤|i|≤d+k‍|Xi|m≤ck,mn(d+k,d).In this proposition, we use the fact mentioned above that whend>k, then the sparse grid X is a union of Xi over i as indicated in the sum given here. Given this result, the following theorem is nearly immediate.Theorem 4.
For fixedk, there is ck>0 so that for d>k, the coefficients in the gPC expansion of A(d+k,d)(f) may be found using (14) with computation time bounded by ckn(d+k,d); that is, for fixed k, the running time is linear in the number of grid points.We have not attempted to calculate the best possibleck in the proof of this result. In fact, numerical results given below show that in some cases, the coefficient ck actually decreases as k increases due to the spread of fixed overhead costs over a small number of points when k and d are small. Additionally, the numerical results show that the time per point evaluated is roughly constant through k=8.Note that this algorithm computes the gPC coefficients of the interpolating polynomial rather than the original function itself. However, for aCr function, f, with r≥1, the error estimate from (5) implies that the interpolating function converges uniformly to f as the depth, k, increases. Since the gPC coefficients are obtained by weighted integration of f against the orthogonal polynomials, the gPC coefficients for the interpolating polynomials converge to the gPC coefficients for f, and (5) provides a means to estimate the error in these coefficients based on the set of orthogonal polynomials and their corresponding weights.Proof of Proposition3.
Given|i|=d+κ between d and d+k, consider the number of points in the corresponding grid Xi. Let r=r(i) be the number of nontrivial entries in i (i.e., entries larger than 1). Suppose without loss that the first r entries of i are nontrivial with entries ij=κj+1≥2. Then, the number of points in this grid is ∏j=1r‍(2κj+1). If κ1 is the only entry larger than 1, then κ1=κ-(r-1), and so the number of points is exactly 3r-1(2κ-r+1+1). A simple counting argument implies that any other distribution of nontrivial entries produces no more grid points, so
(16)|Xi|≤3r-12κ-r+2.
To constructi with |i|=d+κ and r nontrivial entries, we select r out of d nontrivial entries, then distribute κ elements to these entries and so that each entry gets at least one element. There are (dr)(κ-1r-1) ways to do this. Summing from r=1 to κ, and using (dr)≤dr, we have
(17)∑|i|=d+κ‍|Xi|m≤∑r=1κ‍(dr)(κ-1r-1)(3m)r-1(2m)κ-r+2≤(2m)κ+1d∑q=0κ-1‍(κ-1q)(3m)qdq(2m)q.
The summation on the right hand side is a binomial expansion of a power of κ-1, and distributing a factor of (2m)κ-1 into this power converts the right hand side to 4md(3md+2m)κ-1. Using this in place of the sum over i and summing over κ gives
(18)4md∑κ=1k‍(3md+2m)κ-1=4md(3md+2m)k-13md+2m-1≤2dk4m(3m+2md)k-1.
As noted above, n(d+k,d) is asymptotically equal to 2kdk/k! for large d in the sense that for fixed k, limd→∞n(d+k,d)/(2kdk/k!)=1. In particular, there is ak>0 so that dk≤akn(d+k,d) for all d≥1. Hence, the results just given imply
(19)∑κ=1k‍∑|i|=d+κ‍|Xi|m≤2akn(d+k,k)4m(2·3m)k-1.
Hence, taking ck,m=2ak4m(2·3m)k-1, we obtain the desired result.Proof of Theorem4.
Given the sparse grid,X, associated with A(d+k,k) and the vector of values, f, we need first to form the matrices π, D^w, and Φ^ as in (13). From the construction of π, it has exactly one nonzero entry per row and is n×N, where n=n(d,k) is asymptotic to 2kdk/k! as d increases [10], and N=∑d≤|i|≤d+k‍|Xi|. Hence, π and the diagonal matrix, D^w, may be constructed and applied in time O(N).
A given block ofΦ^ is obtained by evaluating products of one-dimensional polynomials on the points of a product grid, Xi. For a nontrivial entry ij=1+κj, the polynomials in xj have degree at most 2κj. Using a three-term recurrence, these polynomials may be evaluated at a given point in O(2κj). Since the projection of Xi to the xj coordinate has 2κj+1 points, these polynomials are evaluated in time O(22κj). This is repeated for each nontrivial κj for a total of O(k22κj), which is O(k|Xi|2). Then for each point in Xi, we multiply at most k polynomials, which are O(k|Xi|). Summing over i and applying Proposition 3, we see that Φ^ may be constructed in time O(kn(d+k,d)).
As noted above, to find the gPC coefficients,c, we define f^=πf, solve Φ^c^=f^ for c^, and then use c=πTD^wc^. Since Φ^ has blocks of size |Xi|×|Xi|, we may solve each block in time at most O(|Xi|3). Combining this with Proposition 3 and the estimates O(N) and O(kN) above, we obtain the theorem.A practical point is that many of the blocks inΦ^ are identical up to a permutation of the rows. Hence, for each of these blocks, we may use a single LU decomposition and appropriate permutations of the entries in f to reduce the total running time (although not, of course, the linear dependence described in the theorem).We note also that the algorithm given here is compatible with the anisotropic adaptive sparse grids of [19]. The analysis given here shows that the time per point depends only on the maximum size of a block in the factorization of Φ-1, and this size is bounded by the maximum depth, k, in both the isotropic and the anisotropic cases.Finally, the analysis for fixedk and increasing d is relevant in that as models become more complex, d increases, while often the important behavior of a model can be captured by a relatively low-order approximation. An example of this is seen in the T-cell signaling model in Section 4.2.


## 4. Numerical results

### 4.1. Efficiency of Conversion

In this section, we provide results of numerical experiments using the algorithm described in the previous sections. All computations were performed in MATLAB 7.7.0 on a Dell Precision PWS690 with an Intel Xeon running at 3 GHz with 3 GB of RAM.To evaluate the running time of the conversion algorithm, we used the test function labeled oscillatory in [10] as
(20)f1(x)=cos(2πw1+∑i=1d‍cixi),
where ci and w1 are chosen at random as indicated in [22]. The domain of definition is [0,1]. The sparse grid for a given dimension and depth was created using the MATLAB package spinterp, version 5.1.1 [23]. The resulting functional values were then used as input to the revised conversion algorithm using the Legendre polynomials as basis.In Figure1, we plot the running time of the conversion algorithm for fixed order of accuracy, k, and increasing dimension. Figure 4(a) clearly shows the linear dependence on the number of points evaluated. Figure 4(b) shows the same data as a function of the number of dimensions. Here, the nonlinear increase for k bigger than 1 is due to the nonlinear dependence of the number of points of evaluation on dimension. Nevertheless, for k=2, the algorithm is reasonably fast even up to 100 dimensions.Time for conversion to Legendre basis as a function of number of points and number of dimensions. Each trace corresponds to a fixed order of accuracy (or depth,k) with the number of dimensions increasing. (a) Time versus number of points showing linear dependence on the number of points evaluated. (b) The same data as a function of number of dimensions (except, k=1 is truncated for scaling reasons). Given the linear relationship on the right, this is essentially a scaled plot of the number of points as a function of the number of dimensions.
(a)
(b)In Figure2, we plot the running time of the conversion algorithm for fixed dimension and increasing order of accuracy, k. Here, the bound in (18) implies that the running time is bounded by c1c2k per point for some constants c1 and c2. However, Figure 4(a) shows that the deviation from nonlinear is relatively small for small k. This is made more precise in Figure 4(b), which shows the time per point as a function of k. Perfect linear dependence would imply that the traces for different dimensions would coincide. While this is essentially true for k≥5 in the data presented, there are significant deviations for small k and d. This is due to fixed overhead time that must average over fewer points in these cases, which gives rise to the decrease in time per point as dimension increases when k≤4. Somewhat more surprisingly, the time per point (after discounting the fixed overhead) is essentially constant up to k=8.Time for conversion to Legendre basis as a function of number of points and order of accuracy (or depth,k). Each trace corresponds to a fixed number of dimensions with the depth increasing. (a) Time versus number of points showing nonlinear dependence on k, particularly in dimensions 3 and 6, for which the maximum depth is 10 and 9, respectively. (b) The same data scaled to give time per point of evaluation and is shown as a function of k. Perfect linear scaling of time with number of points would imply that the time per point is the same for a given value of k, independent of dimension. Apart from fixed overhead costs that skew the results for small k, the time per point is nearly constant over a large range of feasible values of k.
(a)
(b)


### 4.2. Sensitivity Analysis for T-Cell Signaling Model

To illustrate the calculation of global sensitivity coefficients as described in Section2.2, we consider a mathematical model of the T-cell signaling pathway proposed by Lipniacki et al. in [24] and examined further in [25]. This model describes the effect on T-cell activation of agonist and antagonist binding of major histocompatibility complex molecules, with both positive and negative feedback present in the signaling network. We use a form of this model that is a standard ODE system with 19 parameters, 37 state variables, and fixed initial conditions. We focus on the dynamics of phosphorylated ZAP, which plays an important role in positive feedback through the MEK/ERK kinase cascade. We take the uncertain parameter space to be a hypercube with sides given by [pi/5,5pi], where pi is the nominal ith parameter value. For purposes of constructing the gPC expansion and calculating sensitivity coefficients, we transform the parameters to log space. Hence, the function that we evaluate has the form f(t;q), where q is a vector of parameters given by qi=pi5xi, with each xi∈[-1,1], and t∈[0,2000] is time.More precisely, for this example, the nodes of a sparse grid are vectors in[-1,1]19, and each such node gives a corresponding parameter vector q using the transformation just described. For each such q, we evaluate the model at 200 evenly spaced time points between 0 and 2000. For one of these fixed t, say ti, the model evaluations f(ti;q) over all the q given by nodes of a sparse grid allow us to construct an interpolating polynomial, which can then be evaluated at any q in the parameter space, even those that do not arise from nodes of the sparse grid. Of course the resulting interpolation is only an approximation to the true model, but the evaluation of the interpolating polynomial is in this case on the order of 20 times faster than evaluation by integrating the system of ODEs directly, and after change of basis to gPC form, the interpolating polynomial can be used to compute sensitivity coefficients very quickly. Moreover, the remarks in Section 2 quantify the relationship between the number of nodes in the sparse grid and the accuracy of the approximation.Over the selected parameter space, this model has dynamics that can be difficult to capture using interpolation methods. Depending on parameters, the dynamics can exhibit dramatically different time scales, and, moreover, there is a division of parameter space into regions that produce trajectories with relatively high terminal concentrations of pZAP and other regions with relatively low terminal concentrations; very few trajectories have terminal concentrations in an intermediate range. These features can be difficult to capture accurately using polynomial interpolation in the same way that step functions can be difficult to approximate with Fourier series. Figure3 shows a set of sample trajectories of the full model along with trajectories obtained by interpolation using an increasing number of model evaluations. The parameter vectors for these trajectories were chosen randomly from the full parameter space and thus are not sparse grid nodes. These plots show both the difficulty in capturing the dynamics accurately (e.g., the overshoot in the trajectory that rises most quickly) and the fact that interpolation captures the main features of the dynamics even in this difficult case and with accuracy increasing with the number of points.Size sample trajectories for the Lipniacki model, showing the proportion of total ZAP in phosphorylated state (pZAP) versus integration time. (a) shows the computed trajectories. (b)–(d) show interpolated trajectories based on the indicated number of sparse grid nodes.
(a)
(b)
(c)
(d)Sobol’ sensitivity coefficients for each of the 19 parameters as a function of time, based on the indicated number of sparse grid nodes. Although the magnitude of the estimates change slightly between the two plots, the relative ranking is essentially the same.
(a)
(b)After evaluating the function on a sparse grid, we next computed the Sobol’ sensitivity coefficients for each of the 19 parameters as a function of time. As described in Section2.2, this is accomplished by first using the function values f(ti;q) for fixed ti to change to a gPC basis and thereby obtain an interpolating polynomial for f(ti;·); the sensitivity coefficients for this time point are then obtained by summing the squares of coefficients of appropriate terms in this polynomial. This is done for each ti to obtain the sensitivity coefficients as a function of time. The resulting curves are shown in Figure 4, using interpolating polynomials based on 100 function values (a) and 10000 function values (b). As seen there, the estimates of the sensitivity coefficients change in magnitude between the two panels, but the ranking of sensitive versus insensitive is quite consistent between the two. Moreover, this ranking is nearly identical to that given by an independent estimation using a quasi-Monte Carlo method.To demonstrate the utility of this sensitivity analysis, we note that all but 4 of the parameters have sensitivity coefficients that are below 0.1 for the entire time course. The 4 remaining parameters have 3 distinct patterns of influence on the sensitivity. The first pattern is due to the ZAP dephosphorylation rate, which gives rise to a sharp spike in sensitivity near the beginning of the time course (in black in Figure4). This sensitivity at the start of the time course is related to the way in which the model departs from the fixed initial conditions; a larger dephosphorylation leads to a steeper initial decline in the level of pZAP. However, examination of trajectories arising from a range of values for this parameter shows that there is less than 2% variation in these trajectories over the part of the curve corresponding to the peak of sensitivity. Hence, although the sensitivity to this parameter is high in this region, the reason for this sensitivity is not that the trajectories vary greatly in response to changes in this parameter, but rather that changes to the other parameters lead to essentially no change to the output in the early part of the trajectory.The next sensitivity pattern is a moderate hump in the time period of about 50 to about 500 (shown in green and red in Figure4). These are due to the T-cell receptor (TCR) phosphorylation rate (the upper, green curve) and the ZAP phosphorylation rate (the lower, red curve). Not surprisingly, with the remaining parameters held fixed, increased values for either of these rates generally lead to significantly increased levels of pZAP from about time 50 onward.The final sensitivity pattern is shown in Figure4 in the blue curve that rises around time 500, which is due to the SHP phosphorylation rate. SHP plays a role in a negative feedback loop, and this nonlinear dependence gives rise to complex and often long-lasting influences on the trajectories.To indicate the utility of the sensitivity analysis in capturing the dynamics of this model, we took the parameter values used in Figure3, changed all but the 4 parameters discussed as most sensitive back to the nominal values of the original model, and then resimulated. The resulting curves in Figure 5(a) indicate that this simple procedure captures the majority of the dynamical features of the original trajectories of Figure 3(a). We then used sparse grid interpolation on the 4-dimensional grid corresponding to these 4 parameters and produced interpolated curves as in Figure 3. The smaller dimension here means that far fewer points are needed to accurately capture the trajectories relative to the 19 dimensional cases examined earlier.Sample trajectories corresponding to those in Figure3 but with using nominal parameter values for all but the 4 parameters identified as the most sensitive. (a) shows the trajectories obtained by integrating the system of ODEs. (a)–(b) show interpolated trajectories based on the indicated number of points in 4-dimensional sparse grid.
(a)
(b)
(c)
(d)


## 4.1. Efficiency of Conversion

In this section, we provide results of numerical experiments using the algorithm described in the previous sections. All computations were performed in MATLAB 7.7.0 on a Dell Precision PWS690 with an Intel Xeon running at 3 GHz with 3 GB of RAM.To evaluate the running time of the conversion algorithm, we used the test function labeled oscillatory in [10] as
(20)f1(x)=cos(2πw1+∑i=1d‍cixi),
where ci and w1 are chosen at random as indicated in [22]. The domain of definition is [0,1]. The sparse grid for a given dimension and depth was created using the MATLAB package spinterp, version 5.1.1 [23]. The resulting functional values were then used as input to the revised conversion algorithm using the Legendre polynomials as basis.In Figure1, we plot the running time of the conversion algorithm for fixed order of accuracy, k, and increasing dimension. Figure 4(a) clearly shows the linear dependence on the number of points evaluated. Figure 4(b) shows the same data as a function of the number of dimensions. Here, the nonlinear increase for k bigger than 1 is due to the nonlinear dependence of the number of points of evaluation on dimension. Nevertheless, for k=2, the algorithm is reasonably fast even up to 100 dimensions.Time for conversion to Legendre basis as a function of number of points and number of dimensions. Each trace corresponds to a fixed order of accuracy (or depth,k) with the number of dimensions increasing. (a) Time versus number of points showing linear dependence on the number of points evaluated. (b) The same data as a function of number of dimensions (except, k=1 is truncated for scaling reasons). Given the linear relationship on the right, this is essentially a scaled plot of the number of points as a function of the number of dimensions.
(a)
(b)In Figure2, we plot the running time of the conversion algorithm for fixed dimension and increasing order of accuracy, k. Here, the bound in (18) implies that the running time is bounded by c1c2k per point for some constants c1 and c2. However, Figure 4(a) shows that the deviation from nonlinear is relatively small for small k. This is made more precise in Figure 4(b), which shows the time per point as a function of k. Perfect linear dependence would imply that the traces for different dimensions would coincide. While this is essentially true for k≥5 in the data presented, there are significant deviations for small k and d. This is due to fixed overhead time that must average over fewer points in these cases, which gives rise to the decrease in time per point as dimension increases when k≤4. Somewhat more surprisingly, the time per point (after discounting the fixed overhead) is essentially constant up to k=8.Time for conversion to Legendre basis as a function of number of points and order of accuracy (or depth,k). Each trace corresponds to a fixed number of dimensions with the depth increasing. (a) Time versus number of points showing nonlinear dependence on k, particularly in dimensions 3 and 6, for which the maximum depth is 10 and 9, respectively. (b) The same data scaled to give time per point of evaluation and is shown as a function of k. Perfect linear scaling of time with number of points would imply that the time per point is the same for a given value of k, independent of dimension. Apart from fixed overhead costs that skew the results for small k, the time per point is nearly constant over a large range of feasible values of k.
(a)
(b)


## 4.2. Sensitivity Analysis for T-Cell Signaling Model

To illustrate the calculation of global sensitivity coefficients as described in Section2.2, we consider a mathematical model of the T-cell signaling pathway proposed by Lipniacki et al. in [24] and examined further in [25]. This model describes the effect on T-cell activation of agonist and antagonist binding of major histocompatibility complex molecules, with both positive and negative feedback present in the signaling network. We use a form of this model that is a standard ODE system with 19 parameters, 37 state variables, and fixed initial conditions. We focus on the dynamics of phosphorylated ZAP, which plays an important role in positive feedback through the MEK/ERK kinase cascade. We take the uncertain parameter space to be a hypercube with sides given by [pi/5,5pi], where pi is the nominal ith parameter value. For purposes of constructing the gPC expansion and calculating sensitivity coefficients, we transform the parameters to log space. Hence, the function that we evaluate has the form f(t;q), where q is a vector of parameters given by qi=pi5xi, with each xi∈[-1,1], and t∈[0,2000] is time.More precisely, for this example, the nodes of a sparse grid are vectors in[-1,1]19, and each such node gives a corresponding parameter vector q using the transformation just described. For each such q, we evaluate the model at 200 evenly spaced time points between 0 and 2000. For one of these fixed t, say ti, the model evaluations f(ti;q) over all the q given by nodes of a sparse grid allow us to construct an interpolating polynomial, which can then be evaluated at any q in the parameter space, even those that do not arise from nodes of the sparse grid. Of course the resulting interpolation is only an approximation to the true model, but the evaluation of the interpolating polynomial is in this case on the order of 20 times faster than evaluation by integrating the system of ODEs directly, and after change of basis to gPC form, the interpolating polynomial can be used to compute sensitivity coefficients very quickly. Moreover, the remarks in Section 2 quantify the relationship between the number of nodes in the sparse grid and the accuracy of the approximation.Over the selected parameter space, this model has dynamics that can be difficult to capture using interpolation methods. Depending on parameters, the dynamics can exhibit dramatically different time scales, and, moreover, there is a division of parameter space into regions that produce trajectories with relatively high terminal concentrations of pZAP and other regions with relatively low terminal concentrations; very few trajectories have terminal concentrations in an intermediate range. These features can be difficult to capture accurately using polynomial interpolation in the same way that step functions can be difficult to approximate with Fourier series. Figure3 shows a set of sample trajectories of the full model along with trajectories obtained by interpolation using an increasing number of model evaluations. The parameter vectors for these trajectories were chosen randomly from the full parameter space and thus are not sparse grid nodes. These plots show both the difficulty in capturing the dynamics accurately (e.g., the overshoot in the trajectory that rises most quickly) and the fact that interpolation captures the main features of the dynamics even in this difficult case and with accuracy increasing with the number of points.Size sample trajectories for the Lipniacki model, showing the proportion of total ZAP in phosphorylated state (pZAP) versus integration time. (a) shows the computed trajectories. (b)–(d) show interpolated trajectories based on the indicated number of sparse grid nodes.
(a)
(b)
(c)
(d)Sobol’ sensitivity coefficients for each of the 19 parameters as a function of time, based on the indicated number of sparse grid nodes. Although the magnitude of the estimates change slightly between the two plots, the relative ranking is essentially the same.
(a)
(b)After evaluating the function on a sparse grid, we next computed the Sobol’ sensitivity coefficients for each of the 19 parameters as a function of time. As described in Section2.2, this is accomplished by first using the function values f(ti;q) for fixed ti to change to a gPC basis and thereby obtain an interpolating polynomial for f(ti;·); the sensitivity coefficients for this time point are then obtained by summing the squares of coefficients of appropriate terms in this polynomial. This is done for each ti to obtain the sensitivity coefficients as a function of time. The resulting curves are shown in Figure 4, using interpolating polynomials based on 100 function values (a) and 10000 function values (b). As seen there, the estimates of the sensitivity coefficients change in magnitude between the two panels, but the ranking of sensitive versus insensitive is quite consistent between the two. Moreover, this ranking is nearly identical to that given by an independent estimation using a quasi-Monte Carlo method.To demonstrate the utility of this sensitivity analysis, we note that all but 4 of the parameters have sensitivity coefficients that are below 0.1 for the entire time course. The 4 remaining parameters have 3 distinct patterns of influence on the sensitivity. The first pattern is due to the ZAP dephosphorylation rate, which gives rise to a sharp spike in sensitivity near the beginning of the time course (in black in Figure4). This sensitivity at the start of the time course is related to the way in which the model departs from the fixed initial conditions; a larger dephosphorylation leads to a steeper initial decline in the level of pZAP. However, examination of trajectories arising from a range of values for this parameter shows that there is less than 2% variation in these trajectories over the part of the curve corresponding to the peak of sensitivity. Hence, although the sensitivity to this parameter is high in this region, the reason for this sensitivity is not that the trajectories vary greatly in response to changes in this parameter, but rather that changes to the other parameters lead to essentially no change to the output in the early part of the trajectory.The next sensitivity pattern is a moderate hump in the time period of about 50 to about 500 (shown in green and red in Figure4). These are due to the T-cell receptor (TCR) phosphorylation rate (the upper, green curve) and the ZAP phosphorylation rate (the lower, red curve). Not surprisingly, with the remaining parameters held fixed, increased values for either of these rates generally lead to significantly increased levels of pZAP from about time 50 onward.The final sensitivity pattern is shown in Figure4 in the blue curve that rises around time 500, which is due to the SHP phosphorylation rate. SHP plays a role in a negative feedback loop, and this nonlinear dependence gives rise to complex and often long-lasting influences on the trajectories.To indicate the utility of the sensitivity analysis in capturing the dynamics of this model, we took the parameter values used in Figure3, changed all but the 4 parameters discussed as most sensitive back to the nominal values of the original model, and then resimulated. The resulting curves in Figure 5(a) indicate that this simple procedure captures the majority of the dynamical features of the original trajectories of Figure 3(a). We then used sparse grid interpolation on the 4-dimensional grid corresponding to these 4 parameters and produced interpolated curves as in Figure 3. The smaller dimension here means that far fewer points are needed to accurately capture the trajectories relative to the 19 dimensional cases examined earlier.Sample trajectories corresponding to those in Figure3 but with using nominal parameter values for all but the 4 parameters identified as the most sensitive. (a) shows the trajectories obtained by integrating the system of ODEs. (a)–(b) show interpolated trajectories based on the indicated number of points in 4-dimensional sparse grid.
(a)
(b)
(c)
(d)


## 5. Conclusions

The algorithm described here provides an efficient method to convert function values to a gPC approximation to a given function and thereby to estimate global Sobol’ sensitivity coefficients. The method relies on a block-diagonal factorization of the matrix for changing basis from Lagrange polynomials to orthogonal polynomials based on function values at a Smolyak sparse grid. For a fixed degree of accuracy and increasing dimension, this algorithm is linear in the number of points of evaluation. Moreover, for fixed dimension, the time per point of evaluation is nearly constant as the degree of accuracy increases up to aboutk=8.We applied this method to a time-varying model of T-cell signaling and showed that the majority of the dynamics of this model could be explained by variation in only 4 parameters: the ZAP dephosphorylation rate, and the ZAP, TCR, and SHP phosphorylation rates. Taken together, the combination of efficient sensitivity analysis and interpolation provides an easily applied and efficient method for understanding the main dynamics of a complex biological model.



---

*Source: 562767-2013-04-11.xml*