Kriging-Based Design of Experiments for Multi ... - ACS Publications

Research Article pubs.acs.org/journal/ascecg

Kriging-Based Design of Experiments for Multi-Source Exposure− Response Studies in Nanotoxicology Ying Pei,† Feng Yang,*,† Xi Chen,‡ Nianqiang Wu,§ and Kai Wang†,∥

Downloaded via UNIV OF TOLEDO on June 29, 2018 at 20:29:51 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

†

Department of Industrial and Management System Engineering, West Virginia University, P.O. Box 6070, Morgantown, West Virginia 26506, United States ‡ Department of Industrial and Systems Engineering, Virginia Polytechnic Institute and State University, Blackburg, Virginia 24601, United States § Department of Mechanical and Aerospace Engineering, West Virginia University, P.O. Box 6070, Morgantown, West Virginia 26506, United States S Supporting Information *

ABSTRACT: For the sustainable development of nanomaterials (NMs), it is critical to assess the toxic potency of NMs. One of the major challenges with toxicology studies of NMs, compared to traditional materials or chemicals, lies in the large NM variety (or source variations) caused by their various physicochemical properties. How should we efficiently design multi-source biological experiments in nanotoxicolgy, which are extremely time consuming and expensive, to reduce the cost/time for NM toxicity assessment? This work developed a two-stage experimental design procedure based on a statistical model, stochastic kriging with qualitative factors (SKQ). With a given experimental budget, the SKQ-based design method aims at achieving the highest-quality SKQ, which synergistically models the nanotoxicology data from multiple sources (e.g., NM types). The method determines the experimental design (that is, the sampling locations and allocations) in such a way that the resulting sampling data allow SKQ to realize its maximum potential to pool information across multiple sources for efficient modeling. The design procedure also inherits the general advantages of stochastic kriging in the sense that the design is particularly tailored to model the possibly nonlinear and complex relationships and heterogeneous data variances. Through empirical studies, the experimental efficiency of the design procedure for multi-source experiments is demonstrated. KEYWORDS: Risk analysis, Design of experiments, Nanotoxicology, Kriging, Multi-source data

■

INTRODUCTION With the rapid development of nanotechnology, nanomaterials (NMs) find increasing applications in energy,1,2 biomedicine,3−5 and environment.6−9 The production, use, and disposal of NMs inevitably lead to their appearance in air, water, soils, etc., and hence raise safety, health, and environmental concerns. For the safe and sustainable development of nanotechnology, the risk assessment of NMs plays a critical role. One of the most fundamental steps in assessing the risk of an NM (or any other substance) is to characterize the NM by its exposure−response profile,10,11 which describes the dependence of the adverse bioactivity effects (the responses) upon the NM exposure conditions.12 The exposure condition is typically specified through the settings of two factors: NM dosage and time factor involved (e.g., post-exposure time for acute studies). To obtain the characteristic profile of an NM, biological experiments need to be performed at different exposure conditions to observe the corresponding bioactivity responses of animals. Such biological experiments are extremely expensive and time consuming, and the particular challenge with the exposure−response studies of NMs, compared to traditional materials or chemicals, lies in the large NM variety caused by © 2017 American Chemical Society

their various physicochemical properties (e.g., chemical compositions, shape, size, and surface chemistry). How should we design exposure−response experiments across multiple sources (i.e., determine the sampling locations and allocations) for the efficient utilization of limited resources? To address this question, a kriging-based design of experiments (DOE) procedure is developed in this paper. The DOE method is built upon the stochastic kriging with qualitative factors (SKQ) model developed in our previous work.13 Thus, it inherits the advantages of SKQ such as high flexibility and generality and seeks to generate a design that allows SKQ to realize its maximum efficiency in modeling multi-source data. Specifically, the SKQ-based DOE aims at designing multisource experiments to obtain a most informative data ensemble across multiple sources. The data ensemble is most informative in the sense that when being synergistically modeled by SKQ, information will be pooled from all the data sources leading to the highest-quality model (e.g., exposure−response model) for Received: December 8, 2016 Revised: February 25, 2017 Published: February 28, 2017 3223

DOI: 10.1021/acssuschemeng.6b02981 ACS Sustainable Chem. Eng. 2017, 5, 3223−3232

Research Article

ACS Sustainable Chemistry & Engineering

exposure−response profiles of NMs, and its efficiency is illustrated via comparison with the other two DOE methods: the traditional design commonly used in toxicology studies and a space-filling design in the DOE literature.

each source (e.g., NM type). When the target relationships are nonlinear or there is variance heterogeneity in the data, both of which are commonly encountered in nanotoxicology studies, the optimal design that optimizes the model estimation quality is dependent on the true underlying response surfaces and variance patterns, which are unknown. To circumvent this dilemma, the DOE procedure employs a two-stage paradigm: In the first stage, some preliminary experiments are carried out, on which SKQ modeling is performed to derive information regarding the target relationships and data variance; in the second stage, the information obtained from the previous stage is utilized to guide the Stage 2 design, which aims at optimizing the quality of the SKQ models fitted from both stages of data for the target response surfaces. The SKQ-based two-stage design method represents a new addition to the existing literature of experimental design, which can be generally divided into two groups: model-independent versus model-based designs. In these two groups, the design methods most relevant to the current study are briefly reviewed as follows: Both the traditional designs commonly used by biology experimenters and the space-filling designs14,15 fall into the category of model-independent designs, which have no bearing on the models for the target response surfaces. In the current DOE practice for toxicology studies, the adopted traditional designs are typically generated based on empirical experiences in a somewhat arbitrary manner:16,17 They usually involve equally spaced levels (on a linear or log scale) in dose and/or time range, and the same design is typically used across multiple sources. A space-filling design seeks to provide an even coverage of the design regions of interest. In particular, in the presence of multiple sources, sliced space-filling designs18,19 have been developed. In contrast to disregarding the target response surfaces, model-based designs aim at optimizing the quality of the resulting estimated models for the target surfaces. Since our design method is based on SKQ, a kriging model, we focus on reviewing kriging-based designs. Depending on whether or not the response is stochastic, there is deterministic kriging (DK) versus stochastic kriging (SK), based on which respective design methods have been developed. The majority of the work has been on DK-based designs,20−24 which only need to determine the location of design points with one sample assigned to each point. Some research efforts have been devoted to SK-based designs,25,26 which determine not only the design-point locations but also the sample size at each point; SK models the effects of quantitative factors only. As pointed out in Wang et al.,13 SKQ represents the kriging model that models the variability across replications (randomness in responses) and the variability arising from quantitative as well as qualitative factors (e.g., source factors). Accordingly, our SKQ-based design method utilizes the information regarding the target response surfaces and the variance structures to find the design (that is, the design-point locations and sample allocations) for multi-source experiments that leads to the fitted SKQ model of the highest quality. The remainder of this paper is organized as follows: The Statement of the Research Problem section describes in precise terms the problem of designing multi-source experiments for exposure−response studies. A brief review of SKQ and its advantages is given in the Review of Stochastic Kriging with Qualitative Factors (SKQ) section. The SKQ-based two-stage design procedure is detailed in the SKQ-Based Two-Stage Design Procedure section. In the Empirical Studies section, the design procedure is applied to two cases to characterize the

■

STATEMENT OF THE RESEARCH PROBLEM For the exposure−response studies of an NM, biological experiments need to be carried out under a range of experimental conditions. An experimental condition is defined by the combination of a number of factors, which can be divided into two categories, quantitative and qualitative factors. • Quantitative factors typically include but are not limited to the toxicant dosage administered to an animal and the time factor. Depending on the time scope of the toxicology study, the time factor could be exposure time for long-term studies or post-exposure time for acute studies. The vector x is used to represent the quantitative factors considered. • Qualitative factors mainly include the various source factors such as the NM type, the conducting laboratory for experiments, and the exposure routes of animals. The qualitative factors are denoted by the vector z. The experimental condition is specified in terms of the factor vector w = (xT, zT)T. The random response obtained from an animal subject at w can be generally written as @(w) = E[@(w)] + ε(w) = Y(w) + ε(w)

(1)

where Y(w) = E[@(w)] represents the true expected response and ε(w) the random zero-mean error accounting for the variability across replications (or animal subjects). A setting of the qualitative factors z specifies a combination category, say cq, representing one source of data. The total of Q data sources specified by the categories of z are denoted as {cq; q = 1, 2, ..., Q}. The biological data collected at a range of w settings are represented as {(wi , @j(wi)); i = 1, 2, ..., I ; j = 1, 2, ..., n(wi)}

(2)

where I denotes the number of distinct design points (i.e., factor settings at which experiments are performed), wi the ith design point, @j(wi) the response from the jth replication at wi, and n(wi) the number of replications performed at wi. The task of DOE in this work is to determine, with a given experimental budget, the location of design points {wi; i = 1, 2, ..., I} and the sample allocations {n(wi); i = 1, 2, ..., I} for the multi-source experiments so that an SKQ model of the highest quality can be fitted from the sample data approximating the target exposure−response surfaces. As given above, the qualitative vector z may include a range of source factors in nanotoxicology studies. Hence, our method can be used to design experiments involving NMs of different types,27 across different laboratories,28,29 and/or via different animal exposure routes. As mentioned earlier, the DOE is built in a two-stage learning framework to enable the search of the optimal design for the unknown target response surfaces. Such a two-stage framework is also consistent with the practice of biologists: In biological studies, it is typical to perform some preliminary experiments before launching the relatively large-scale experimentation, and the former are usually used to assist the design of the latter. The two-stage framework can be straightforwardly extended to include multiple stages: All the existing relevant 3224


Research Article

ACS Sustainable Chemistry & Engineering data (including those in the literature) can be considered as the Stage 1 data and used to guide the design of a new batch of experiments when a follow-up stage of experimentation is to be carried out.

For I distinct design points, the I × 1 vector of sample averages is denoted as

REVIEW OF STOCHASTIC KRIGING WITH QUALITATIVE FACTORS (SKQ) The DOE method in this work is developed based on the SKQ modeling of multi-source data. Hence, a brief review of the SKQ developed in Wang et al.13 is provided herein. As pointed out in Wang et al.,13 SKQ is the first kriging model that is able to accommodate the variability arising from quantitative as well as qualitative factors and the variability across replications. SKQ models the dependence of a continuous response upon the factors w = (xT, zT)T, with x = (x1 , x 2 , ..., xd)T ∈ d and z = (z1, z2, ..., zL)T. There are L qualitative factors, and each factor zS(S = 1, 2, ..., L) has a number of category levels. The response at w from the jth replication (animal subject) is modeled by SKQ as

The vector of sample average errors is represented as

Y ̅ = (Y ̅ (w1), Y ̅ (w2), ..., Y ̅ (wI ))T

■

@j(w) = Y(w) + εj(w) = f(w)T β + M(w) + εj(w)

ε ̅ = (ε ̅ (w1), ε ̅ (w2), ..., ε ̅ (wI ))T

Cov[M(w), M(w′)] = δ 2 × Corr[M(w), M(w′)] L

= δ × [∏ τz(SS,)zl′] × K (x , x′) 2

S= 1

n(wi)

∑ @j(wi) = β0 + M(wi) + j=1

1 n(wi)

(6)

where δ2 is the variance of the Gaussian process. The correlation Corr[M(w), M(w′)] is decomposed into two parts: ΠSL= 1τz(SS,)zl′ and K(x, x′). For the estimation of an SKQ, specific functional forms need to be assumed for both parts. The correlation across the quantitative settings is represented by K(x, x′), for which a range of functional forms are provided in the literature.30,31The most widely adopted function is the exponential correlation function

(3)

The expectation Y(w) consists of two parts: f(w) β and M(w). Here, f(w) is a vector of known functions of w and β a vector of unknown parameters of compatible dimension. In this work, we set f(w)Tβ = β0, which has been widely accepted as sufficient for most applications. The term M(w) represents a mean-zero stationary Gaussian process and intends to capture the variability due to the factors w, which is referred to as the extrinsic variability. The intrinsic variability refers to the randomness of ε(w). The random noise {εj(w); j = 1, 2, ...} at w is assumed to have mean zero and be independent and identically distributed (i.i.d.) across replications. The error variance Var[ε(w)] can be w dependent. With the sample data (eq 2), the sample average of the responses at wi across the n(wi) replications is obtained as 1 n(wi)

(5)

n(wi) −1 with ε(w ̅ i) = n(wi) ∑j = 1 εj(wi), i = 1,2, ..., I. The extrinsic variability is modeled in SKQ by M(w), which is specified by its covariance

T

Y ̅ (wi) =

(4)

d

K (x , x′) = exp{∑ −θh|xh − xh′ | p } h=1

(7)

where θ = (θ1, θ2, ..., θd) is a vector of unknown parameters. In eq 6, the term ΠSL= 1τz(SS,)zl′ models the correlations across different categories of qualitative factors, and the vector Φ denotes the unknown parameters involved in the cross-category correlation model. Choices of the functional forms for τz(SS,)zl′ were discussed in details in Wang et al.13 Given the data (eq 2) collected at I distinct design points, the I × I variance−covariance matrix ΣM is defined as

n(wi)

∑ εj(wi) j=1

Σ M = δ 2 × R(θ , Φ)= ⎛ 1 Corr[M(w1), M(w2)] ⎜ ⎜ Corr[M( w ), M( w )] 1 2 1 δ2 × ⎜ ⋮ ⋮ ⎜ ⎜ ⎝ Corr[M(wI ), M(w1)] Corr[M(wI ), M(w2)]

where R(θ,Φ) denotes the correlation matrix with each element representing a correlation. Each element correlation can be decomposed into two parts as explained above and involves the unknown parameters θ and Φ. For an arbitrary w0, the I × 1 vector ΣM(w0, ·) is defined as

⋯ Corr[M(w1), M(wI )]⎞ ⎟ ⋯ Corr[M(w2), M(wI )]⎟ ⎟ ⋱ ⋮ ⎟ ⎟ ⋯ 1 ⎠

(8)

with v(w0, θ, Φ) being a correlation vector involving w0, θ, and Φ. Denote Σε as the I × I variance−covariance matrix of vector ε defined in eq 5. With i.i.d. random errors, Σε is a I × I diagonal matrix

⎛ Corr[M(w0), M(w1)]⎞ ⎟ ⎜ ⎟ ⎜ Corr[M( ), M( )] w w 0 2 Σ M(w0, ·) = δ 2 v(w0, θ , Φ) = δ 2⎜ ⎟ ⋮ ⎟ ⎜ ⎟ ⎜ ⎝ Corr[M(w0), M(wI )]⎠

Σε = diag{Var[ε(w1)]/n(w1), Var[ε(w2)]/n(w2), ..., Var[ε(wI )]/n(wI )}

(10)

For a data set (eq 2), the SKQ estimation and inference procedure in Wang et al.13 are summarized as follows: (1) Estimate Σε as

(9) 3225


Research Article


(exposure−response surfaces for multiple sources) estimated from both stages of data. Initial Design and Preliminary Modeling. At Stage 1, N1 samples are to be assigned. Denote the initial design as

Σ̂ε = diag{Var[̂ ε(w1)]/n(w1), Var[̂ ε(w2)]/n(w2), ..., Var[̂ ε(wI )]/n(wI )} (11)

where Var[̂ ε(wi)] =

1 n(wi) − 1

⎛ w1 w2 ⋯ wI1 ⎞ ⎟ D1 = ⎜ ⎝ n1 n2 ⋯ n I1 ⎠

n(wi)

∑ (@j(wi) − Y ̅(wi))2 , i = 1, 2, j=1

..., I

where I1 is the number of distinct design points. Given the sample size N1, the initial design (eq 15) is determined as follows: At this initial stage, ni is set as n for i = 1, 2, ..., I1, with absent information on the variance pattern. The value of n is specified by the experimenter based on the experience with biological experiments with NMs, which is typically done to determine a traditional design in the current practice. The number of distinct design points I1 is then calculated as I1 = N1/n. The I1 sampling points will then be allocated to the Q slices (sources) of the region χ. The purpose of the initial design is to provide a fair coverage of the design space to gain some information on both the target surfaces and the variance structure throughout the multi-slice design regions, and a range of designs can be adopted for this initial stage. It could be the traditional design that biologists have been using. In each dose−time region for one source, the same evenly spaced design points are selected for sampling; control points are typically included in such designs. In this work, Sliced Latin Hypercube Design (SLHD)18 is adopted for the Stage 1 design as an example. SLHD is a special space-filling design that can be partitioned into slices of smaller Latin hypercube designs (LHD). Each slice of the design achieves maximum uniformity in any one-dimensional projection and provides an even coverage of χ for each slice. When collapsed across the slices, all the sampling points in χ have the maximum stratification in any one-dimensional projection. Although no controls were included in the SLHD for the empirical cases (Empirical Studies section), SLHD can be easily adapted to include controls if needed by rounding the lowest-dose point to be a control point with zero dose. From the initial data, which are denoted as the following, {(wi , @j(wi)); i = 1, 2, ..., I1; j = 1, 2, ..., n(wi)}, two preliminary models are estimated. First, the SKQ estimation (Review of Stochastic Kriging with Qualitative Factors (SKQ) section ) is performed quantifying the relationship between the expected response Y(w) and w (that is, the multiple response surfaces across sources). Second, deterministic kriging34 is carried out based on the data pair {(wi , Var[̂ ε(wi)]); i = 1, 2, ..., I1}, where Var[̂ ε(wi)] represents the sample variances as calculated in eq 12. The second kriging model approximates Var[ε(w)] as a function of w. These two kriging models, which contain the target surface and variance information derived from the initial data, are utilized to guide the follow-up design in Stage 2. Design Augmentation. At Stage 2, the task is to determine the design for the N2 = N − N1 samples, which are to be assigned to I2 = I − I1 distinct design points. The Stage 2 design is denoted as

(12)

(2) Replace Σε by Σ̂ε in the maximum likelihood (ML) function and obtain the ML estimates (β̂0, δ̂2, θ̂, Φ̂ ), which fully specify the SKQ model. (3) For an arbitrary w0, the expected response Y(w0) can be estimated by 2 Y(̂ w0) = β0̂ + v(w0, θ̂ , Φ̂ )T [δ ̂ R(θ̂ , Φ̂ ) + Σ̂ε]−1

(Y ̅ − β0̂ 1I )

(13)

The mean squared error (MSE) of Ŷ (w0) is obtained as ̂ ̂ w0)] = δ 2̂ − δ 4̂ v(w0, θ ̂, Φ̂)T [δ 2̂ R(θ ̂, Φ̂) + Σ̂ε]−1 MSE[Y( v(w0, θ ̂, Φ̂) 2 + η2(1TI [δ ̂ R(θ ̂, Φ̂) + Σ̂ε]−1 1I )−1

(14)

■

where η = 1 −

1TI

(15)

[δ R(θ̂, Φ̂ )+ Σ̂ε]−1 v(w0, θ̂, Φ̂ )δ2̂ . 2̂

SKQ-BASED TWO-STAGE DESIGN PROCEDURE An SKQ-based two-stage procedure is developed to determine the design {wi, n(wi); i = 1, 2, ..., I} for multi-source experiments. The prespecified inputs required by the procedure are given as follows: • N: total number of samples (animals) to be allocated to all the sources considered. This number depends on the experimental budget available. • N1: number of samples assigned to Stage 1, which implies that the Stage 2 sample size is N − N1. It is recommended to set N1 as 1/4−1/2 of N.32 • H: feasible space of the variable vector w. Here, H includes Q slices (Q sources) of the region χ, denoting the feasible region of the quantitative vector x. • I: total number of distinct design points to be assigned. In biological experiments, it is typical to assign at least three samples to a design point, and thus, I ≤ N/3. Also, I is recommended to be set at least about 10 times the dimension of χ33 to provide an adequate coverage of the region χ. An overview of the SKQ-based two-stage procedure is briefed as follows: In Stage 1, N1 pilot experiments are performed following the initial design, which is discussed in the Initial Design and Preliminary Modeling section. On the basis of the preliminary data, two kriging models will be fitted: (i) one quantifies the dependence of the expected response Y(w) upon w, and (ii) the other approximates the variance Var[@(w)] as a function of w. In Stage 2, the information contained in the two fitted models is utilized to design the remaining N − N1 experiments. The follow-up design aims at optimizing the quality of the resulting response surfaces

⎛ wI1+ 1 wI1+ 2 ⋯ wI ⎞ ⎟ D2 = ⎜ ⎝ n I1+ 1 n I1+ 2 ⋯ nI ⎠

(16)

and is determined in such a way that the quality of the final SKQ, which is fitted from both stages of data to model the multi-source response surfaces, is optimized. 3226


Research Article


Figure 1. Two-stage experimental design: stars denote the Stage 1 design points, and dots denote the Stage 2 design points.

Hence, for a candidate design of D2, the corresponding IMSE criterion can be evaluated based on the prior obtained estimates for (β0, δ2, θ, Φ), the variance model, and D1 determined in Stage 1. This evaluation ability provides the necessary basis to perform a numeric search for the optimal values of D2. Solving eq 17 is challenging due to the high dimension of the decision variables involved in D2 and the integer nature of ni (i = I1 + 1, I1 + 2, ..., I). We approach this problem by utilizing the approximate relationships between the optimal values of {n(wi); i = 1, 2, ...} and the given design locations {wi; i = 1, 2, ...} and by an iterative procedure. Given the remaining sample size N2, the first stage design D1 and the design-point locations of the second stage {wi; i = I1 + 1, I1 + 2, ..., I}, the value of {n(wi); i = I1 + 1, I1 + 2, ..., I} that minimizes the IMSE of the SKQ model, which is fitted from the N-sample data following the design {(wi, n(wi)); i = 1, 2, ..., I}, can be approximated as

The quality of a model can be measured by a range of metrics such as the integrated mean squared error (IMSE) and the comparison index based criterion.26 An appropriate performance measure can be selected depending on the modeling purpose. Herein, as an example, IMSE is adopted as the design criterion, and the design optimization problem in Stage 2 can be formulated as MinimizeD2 IMSE =

∫w ∈H MSE[Y ̂(w0)]d w0 0

(17)

I

Subject to

∑

ni = N2

i = I1+ 1

(18)

where I1, I, and N2 are all predetermined parameters. The objective criterion IMSE represents the integration of MSE[Ŷ (w0)], and the MSE of the response estimates Ŷ (w0) over the feasible space H. The estimate Ŷ (w0) at an arbitrary w0 is provided by the SKQ model fitted from both stages of data following the design D = D1 ∪ D2. How is the criterion IMSE related to the design D? Recall that for a fitted SKQ model, the MSE of the response estimates MSE[Ŷ (w0)] can be estimated as eq 14, which is rewritten as follows for convenience of discussions:

n(wi) ≈ N2

VC i i I ∑j=1

VjCj

(20)

The detailed derivation is given in the Replication versus Location Dependence section of the Supporting Information. −1 In eq 20, Vi = Var[ε(wi)] and Ci = [Σ−1 M WΣM ]ii, where W is the I × I matrix with elements Wij = ∫ w0∈Hr0ir0jdw0 and r0i = Corr[M(wi), M(w0)]. For a candidate design, both Vi and Ci can be estimated based on the fitted models from Stage 1. Recognizing eq 20, an iterative procedure is developed as follows to solve the design optimization (eq 17): Inputs: (a) D1, the Stage 1 design (b) SKQ model fitted from the Stage 1 data approximating the target multi-source response surfaces (c) Kriging model fitted from the Stage 1 data approximating the random error variance Var[ε(w)] as a function of w

2 4 2 MSE[Ŷ ̂ (w0)] = δ ̂ − δ ̂ v(w0, θ̂ , Φ̂ )T [δ ̂ R(θ̂ , Φ̂ ) + Σ̂ε]−1 v(w0, θ̂ , Φ̂ ) 2 + η2(1TI [δ ̂ R(θ̂ , Φ̂ ) + Σ̂ε]−1 1I )−1 ,

(19)

where η = 1 − 1TI [δ2̂ R(θ̂, Φ̂ ) + Σ̂ε]−1 v(w0, θ̂, Φ̂ )δ̂2. As shown from eq 19, an MSE (and hence the IMSE) depends on three items: the SKQ parameters (β0, δ2, θ, Φ) for the target response surfaces, which have been estimated from Stage 1 and can be plugged into eq 19; the variance model specifying the relationship between Var[ε(w)] and w, which have been estimated from Stage 1 as well; and the design D = D1 ∪ D2, with D1 given from Stage 1 and D2 to be determined via the optimization. More specifically, in eq 19, both R(θ̂, Φ̂ ) and v(w0, θ̂, Φ̂ ), which are defined in eqs 8 and 9, respectively, depend on the design D and the SKQ parameters (β0, δ2, θ, Φ). The diagonal variance matrix Σ̂ε, which is defined in eq 10, has each element as Var[ε(wi)]/n(wi), i = 1, 2, ..., I and is dependent on D; the functional dependence of Var[ε(w)] upon w has been approximately established by the variance model.

Step 0: Set ni ≈ N2

VC i i I

∑ j = 1 VjCj

for i = I1 + 1, I1 + 2, ..., I; solve

eq 17, with respect to {wi; i = I1 + 1, I1 + 2, ..., I}. Step 1: Fix {wi; i = i = I1 + 1, I1 + 2, ..., I} at its most recently obtained values, and solve eq 17 with respect to {n(wi); i = I1 + 1, I1 + 2, ..., I}. Step 2: Fix {n(wi); i = I1 + 1, I1 + 2, ..., I} at its most recently obtained values, and solve eq 17 with respect to {wi; i = I1 + 1, I1 + 2, ..., I}. 3227


Research Article


factor vector w = (xT, z)T, and the response of interest is BAL (bronchoalveolar lavage) PMNs measured in the units of 103/ mouse. Case 1: Simulation Model. A Monte Carlo simulation model is developed to generate dose−time−response data mimicking the real experimental data in the TiO2 toxicology study described above. The simulation model is specified as follows: The true expected dose−time−responses for the two sources (short and long nanobelts) are represented as {Y(x, c1), Y(x, c2)}, with specific expressions given by models S12 and S13 in the The Simulation Model of Case 1 section of the Supporting Information. The error variance is dependent on w, and the true variance models are given as

Step 3: Repeat Steps 1 and 2 until there is no significant changes in {wi; i = I1 + 1, I1 + 2, ..., I} or {n(wi); i = I1 + 1, I1 + 2, ..., I}. Through the iterative procedure, the size of the optimization problems to be solved is reduced by half. In our empirical experience, it suffices to perform one round of Steps 1 and 2 to achieve convergence in D2. The Particle Swarm Optimization (PSO) is employed to solve the optimization problems involved in this procedure. Compared to other global optimization algorithms such as Genetic Algorithm (GA), the reasons of adopting PSO are as follows: First, PSO has significantly better computational efficiency than GA for the same effectiveness35 since it has no evolution operators (crossover and mutation) while GA does. Second, PSO suits high-dimensional problems, which are frequently encountered in the design of multi-source experiments since the momentum effects on particle movement in PSO allows for good diversity and explorations.

■

EMPIRICAL STUDIES The two-stage design procedure is evaluated through two simulation cases, which mimic real exposure−response studies in nanotoxicology.

Var[ε(x, c1)] = (0.33Y(x , c1)0.26 )2

(21)

Var[ε(x , c 2)] = (0.77 exp(Y(x, c 2) × 0.0076))2

(22)

For a source cq (q = 1, 2) and at an exposure level x0, a random response y0 is simulated as y0 = Y(x 0, cq) +

Var[ε(x 0, cq)] × ϵ; q = 1, 2

(23)

where ϵ is a random variable provided by a standard normal random generator.36 The simulation model serves two purposes in this study. First, it is used to generate dose−time−response data via computer experiments for modeling. Second, the true expected response surfaces {Y(x, c1), Y(x, c2)}, which are part of the simulation model, provide the true benchmark to evaluate the SKQ model fitted from data following a certain design and hence to compare different design methods in terms of their efficiency. Case 1: Applying the Two-Stage Procedure. By using the simulation model (Case 1: Simulation Model section) as the sampling approach for data generation, we applied the SKQbased procedure to design and model the multi-source experiments for the TiO2 toxicity study. The inputs of the procedure (listed at the beginning of the SKQ-Based Two-Stage Design Procedure section) are given as follows for this case: N = 240; N1 = 120; H includes Q = 2 slices of the dose−time region χ = [0, 15] × [1, 112]; I = 24. Stage 1: A total of N1 = 120 samples are assigned to the initial stage. The number of replications at each design point is selected to be 10, and thus, the number of distinct design points is I1 = 12. The locations of these 12 points are determined by applying SLHD to the two slices of region χ. A realization of the SLHD space-filling design is plotted as stars in Figure 1(a) and (b), with the number next to each star representing the number of replications assigned to that design point. Each of the two sources (short and long nanobelts) has been assigned 6 points. Thus, the Stage 1 design D1 consists of 12 points (stars in Figure 1) with 10 replications at each point. Following the initial design, model 23 is used to generate 10 i.i.d. random responses at each distinct design point. The initial data set is denoted as {(wi , @j(w)); i = 1, 2, ..., 12; j = 1, 2, ..., 10}, on which normalization is performed so that both the quantitative factors and responses fall within the range of [0, 1]. On the basis of the initial data, two models are fitted: The fitted SKQ model Ŷ (w) approximating the two target dose−time− response surfaces is obtained by applying the SKQ estimation

Figure 2. Traditional design: design points in the dose−time region for each source (short or long nanobelts).

Figure 3. Space-filling design: triangles denote the design points for short nanobelts, and squares denote the design points for long nanobelts.

Case 1. This case is derived from the in vivo study of TiO2 NMs performed in Porter et al.27 The experimental condition is defined by two quantitative factors x = (x1, x2) and one qualitative factor z. The TiO2 dosage administered to an animal is denoted by x1, with x1 ∈ [0, 15] μg, and the post-exposure time is denoted by x2, with x2 ∈ [1, 112] days. The factor z has two categories (sources) {c1, c2}, corresponding to the two different shapes of NMs: c1 for short and c2 for long TiO2 nanobelts. Thus, the experimental condition is specified by the 3228


Research Article


Figure 4. ERMSE box plots for each of the three design methods in Case 1.

Table 1. Twelve Sources {ci; i = 1, 2, ..., 12} Specified by Three Qualitative Factors source

z1

z2

z3

c1 c2 c3 c4 c5 c6 c7 c8 c9 c10 c11 c12

mice rats mice rats mice rats mice rats mice rats mice rats

crystal crystal noncrystal noncrystal crystal crystal noncrystal noncrystal crystal crystal noncrystal noncrystal

belt belt belt belt particle particle particle particle tube tube tube tube

Figure 1 represent one possible SLHD space-filling design for Stage 1. The Stage 2 design is dependent on both the design and randomly sampled data in Stage 1. The final SKQ model is fitted from the randomly sampled data following the two-stage design. Case 1: Evaluation and Comparison. In light of the stochastic nature of its outcome, the two-stage design procedure is statistically evaluated and compared to the other two design methods based on large macroreplications. The two alternative design methods considered are described as follows: (i) The traditional design includes 12 evenly spaced sampling points (Figure 2) in the dose−time region for each of the two sources. Ten replications are assigned to each design point, with a total sample size of 10 × 12 × 2 = 240. (ii) The SLHD space-filling design is performed with 24 distinct design points assigned to the two slices of dose−time region, and the total sample size is also 240 with 10 replications at each point. An example of the generated SLHD design is given in Figure 3, with the triangles representing the design points for the short nanobelts and the squares those for the long nanobelts. Clearly, the three design methods have the same experimental budget of 240 samples to two sources. For each design method, M = 200 macroreplications (that is, independent applications of the method) have been performed as follows: Step 0: Initialize the index m = 1. Step 1: Apply the design method and follow the design to carry out the experiments. On the basis of the collected data, fit the SKQ model approximating the exposure−response surfaces for both short and long nanobelts; denote the fitted SKQ as Ŷ (m)(w). Step 2: Evaluate the goodness of Ŷ (m)(w) by calculating the estimated root mean squared error (ERMSE)

procedure (Review of Stochastic Kriging with Qualitative Factors (SKQ) section). The variance model Var(̂ w) is fitted from the following sample variance data via deterministic kriging: {(wi , Var[̂ ε(w)]); i = 1, 2, ..., 12}. Stage 2: Here, N2 = N − N1 = 120 samples are to be assigned to I2 = I − I1 = 12 distinct design points. The location of these 12 design points and the sample allocation are obtained by solving the design optimization problem (eq 17) using the iterative procedure in the Design Augmentation section. The inputs needed for the iterative procedure have been obtained from Stage 1: Stage 1 design D1, SKQ model Ŷ (w), and variance model Var[̂ w] that are both fitted from the Stage 1 data. The design points obtained from the design optimization for Stage 2 are plotted as dots in Figure 1(a) and (b), with the number next to each dot representing the sample size assigned to that design point. Among the 12 design points, five are assigned to the short nanobelts and seven to the long nanobelts. The samples are not evenly distributed to the 12 design points: the points with higher error variance tend to have more replications assigned to them. Following the Stage 2 design, a new batch of experiments are carried out. From the two stages of data, SKQ is performed to simultaneously model the exposure−response surfaces for both short and long TiO2 nanobelts. It is worth noting that the outcome (the resulting two-stage design, sample data, and fitted SKQ from both stages of data) of applying the two-stage procedure is random: The stars in

ERMSE(m)(*cq) =

1 #[*cq]

(Y ̂

∑

(m)

(w) − Y (m)(w))2 ; q = 1, 2

w ∈ *cq

(24)

In eq 24, *cq represents the collection of 8456 check points in the dose−time region for short (q = 1) or long (q = 2) nanobelts. These check points are evenly spaced dense grids in the dose−time region. The total number of check points in the set *cq is denoted as #[*cq] and is equal to 8456 for both q values. ERMSE measures the average deviation of Ŷ (·) from its true correspondence Y(·) at the check points in * . The true 3229


Research Article


Figure 5. ERMSE box plots for each of the three design methods in Case 2.

includes the two box plots for each of the three design methods. In comparison, the boxes of our two-stage procedure are the lowest, illustrating its best performance in terms of ERMSE among the three methods; the boxes are relatively narrow, showing that the two-stage procedure leads to designs with consistently low ERMSE. The ERMSE performance of the space-filling design is quite unstable, with the two widest boxes. The ERMSE of the traditional design is consistently high, corresponding to the highest and narrowest boxes in the figure. Comparing the boxes of the traditional design and our two-

value Y(·) is available from the simulation models S12 and S13 in the SI. Step 3: If m < M, then m = m + 1, and repeat Steps 1 and 2. Otherwise, terminate the procedure. For each design method, the following estimated root mean squared errors, {ERMSE(m)(*c1); m = 1, 2, ..., 200} and {ERMSE(m)(*c2); m = 1, 2, ..., 200}, are obtained from 200 macroreplications and used to generate the two box plots for short (c1) and long (c2) nanobelts, respectively. Figure 4 3230


Research Article


exposure−response studies of NMs. The design method is closely coupled with stochastic kriging with qualitative factors (SKQ), a statistical model which synergistically models multisource data by pooling information across sources. On the basis of SKQ, the DOE procedure seeks to find the experimental design (the sampling locations and allocations) that enables SKQ to maximize its information-pooling capability and thus to render fitted relationships (e.g., exposure−response) of the highest quality. The DOE is built in a two-stage framework to allow for the learning of the possibly nonlinear relationships and heterogeneous variance structures, and the learned information is utilized to guide the optimal experimental design. The SKQ-based two-stage procedure has been applied to two simulation cases mimicking the exposure−response study of NMs, and its efficiency over the two alternative design methods has been demonstrated. The DOE method developed in this paper is expected to substantially reduce the cost and time in exposure−response studies for nanotoxicology, accelerate the progress toward quantifying the risk, safety, and health effects of NM exposure, support the sustainable development of NMs, and facilitate the safe commercialization of NMs and their products.

stage procedure, it is evident that the latter is about four times more efficient than the former. Case 2. This case is derived from exposure−response studies of Multiwall Carbon Nanotubes (MWCNT) with different properties and performed on different animals. The experimental condition w = (xT, zT)T is specified by two quantitative factors x = (x1, x2) and three qualitative factors z = (z1, z2, z3). Here, x1 denotes the lung-weight normalized dosage deposited in a rodent, x1 ∈ [0,30] μg/g, and x2 is the postexposure time with x2 ∈ [0, 40] days. Also, z1 represents the rodent species including two categories: mice or rats, and z2 is a binary factor denoting whether the NMs are of crystal structure: YES or NO. Here, z3 has three categories corresponding to the three different shapes of NMs: belt, particle, and tube. The response of interest is the PMNs in the cell samples. Case 2: Simulation Model. In this case, there are 12 sources, with each corresponding to a combination formed by the three qualitative factors (Table 1). As in Case 1, the true expected dose−time−responses and variance models for each of the 12 sources are denoted as {Y(x, ci), i = 1, 2, ..., 12} and {Var[ε(x, ci)], i = 1, 2, ..., 12}. The details of the simulation model is given in the The Simulation Model of Case 2 section of the Supporting Information. Case 2: Results and Comparison. The two-stage design procedure was applied to this case and compared with the two alternative designs. The inputs for the two-stage design is provided as follows: The total number of samples available to the 12 sources is N = 864. The number of samples assigned to the initial stage is N1 = 384, and the number of replications in Stage 1 is set as 8 for each design point. Here, H includes Q = 12 slices of the dose−time regions, and the total number of distinct design points I = 108. Under these settings, the design procedure generates a design including I = 108 distinct design points in the Q = 12 dose−time regions. The two alternative designs are specified as follows given the same total sample size N = 864: (i) The traditional design includes nine points in each of the 12 dose−time regions: three evenly spaced levels in both dose and time ranges lead to 3 × 3 = 9 points. At each point, eight replications are assigned, and hence, N = 8 × 9 × 12 = 864. The location of design points remains the same across different macroreplications. (ii) The space-filling design assigns nine design points to each of the 12 slices, with eight replications at each point. Thus, its total samples N = 8 × 9 × 12 = 864 as well. The outcome of this design is random and varies across macroreplications. For each of these three design methods, the following occurs: M = 100 macroreplications have been performed, and {ERMSE(m)(*ci); i = 1, 2, ..., 12; m = 1, 2, ..., 100} are gathered for each design method to generate the 12 box plots for the 12 sources (Figure 5). Compared to the other two design methods, the ERMSE boxes of the two-stage design procedure are the lowest: Given the same experimental budget, the resulting design from the two-stage procedure leads to fitted response surface models of the lowest ERMSE (highest quality). Also, the boxes of the two-stage design are the narrowest, showing its stable performance across different macroreplications.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acssuschemeng.6b02981. Approximate relationship between design-point locations and sample allocations and simulation models of Case 1 and Case 2. (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: 304-293-9477. Fax: 304-293-4970. ORCID

Feng Yang: 0000-0002-6917-7053 Nianqiang Wu: 0000-0002-8888-2444 Present Address ∥

(K.W.): Fifth Third Bancorp 5001 Kingsley Dr, Cincinnati, OH, 45227. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research was partially supported by the National Science Foundation (NSF) Grant CBET-1065931 and National Institute for Occupational Safety and Health (NIOSH). The findings and conclusions in this report are those of the authors and do not necessarily represent the views of NSF or NIOSH. The authors gratefully thank Dr. Eileen Kuempel from NIOSH for valuable discussions.

■

■

REFERENCES

(1) Zhi, M.; Manivannan, A.; Meng, F.; Wu, N. Highly conductive electrospun carbon nanofiber/MnO 2 coaxial nano-cables for high energy and power density supercapacitors. J. Power Sources 2012, 208, 345−353.

SUMMARY A two-stage design of experiments (DOE) procedure was developed to efficiently design multisource experiments in 3231


Research Article

ACS Sustainable Chemistry & Engineering (2) Dai, L.; Chang, D. W.; Baek, J.-B.; Lu, W. Carbon nanomaterials for advanced energy conversion and storage. Small 2012, 8, 1130− 1166. (3) Li, M.; Zhang, J.; Suri, S.; Sooter, L. J.; Ma, D.; Wu, N. Detection of adenosine triphosphate with an aptamer biosensor based on surfaceenhanced Raman scattering. Anal. Chem. 2012, 84, 2837−2842. (4) Ren, N.; Li, R.; Chen, L.; Wang, G.; Liu, D.; Wang, Y.; Zheng, L.; Tang, W.; Yu, X.; Jiang, H.; Liu, H.; Wu, N. In situ construction of a titanate-silver nanoparticle-titanate sandwich nanostructure on a metallic titanium surface for bacteriostatic and biocompatible implants. J. Mater. Chem. 2012, 22, 19151−19160. (5) Giersig, M.; Khomutov, G. B. Nanomaterials for Application in Medicine and Biology; Springer, 2008. (6) Li, M.; Wang, Q.; Shi, X.; Hornak, L. A.; Wu, N. Detection of mercury (II) by quantum dot/DNA/gold nanoparticle ensemble based nanosensor via nanometal surface energy transfer. Anal. Chem. 2011, 83, 7061−7065. (7) Wu, N.; Zhao, M.; Zheng, J.-G.; Jiang, C.; Myers, B.; Li, S.; Chyu, M.; Mao, S. X. Porous CuO-ZnO nanocomposite for sensing electrode of high-temperature CO solid-state electrochemical sensor. Nanotechnology 2005, 16, 2878. (8) Wu, N.; Chen, Z.; Xu, J.; Chyu, M.; Mao, S. X. Impedance-metric Pt/YSZ/Au-Ga 2 O 3 sensor for CO detection at high temperature. Sens. Actuators, B 2005, 110, 49−53. (9) Farré, M.; Gajda-Schrantz, K.; Kantiani, L.; Barceló, D. Ecotoxicity and analysis of nanomaterials in the aquatic environment. Anal. Bioanal. Chem. 2009, 393, 81−95. (10) Edler, L.; Poirier, K.; Dourson, M.; Kleiner, J.; Mileson, B.; Nordmann, H.; Renwick, A.; Slob, W.; Walton, K.; Würtzen, G. Mathematical modelling and quantitative methods. Food Chem. Toxicol. 2002, 40, 283−326. (11) Bretz, F.; Hsu, J.; Pinheiro, J.; Liu, Y. Dose finding-a challenge in statistics. Biom. J. 2008, 50, 480−504. (12) Benchmark Dose Technical Guidance; U.S. Environmental Protection Agency, Washington, DC, 2012. (13) Wang, K.; Chen, X.; Yang, F.; Porter, D. W.; Wu, N. A new stochastic kriging method for modeling multi-source exposureresponse data in toxicology studies. ACS Sustainable Chem. Eng. 2014, 2, 1581−1591. (14) Qian, P. Z.; Tang, B.; Jeff Wu, C. Nested space-filling designs for computer experiments with two levels of accuracy. Statistica Sinica 2009, 19, 287. (15) Qian, P. Z.; Ai, M.; Wu, C. J. Construction of nested spacefilling designs. Annals of Statistics 2009, 37, 3616−3643. (16) Hamilton, R.; Wu, N.; Porter, D.; Buford, M.; Wolfarth, M.; Holian, A. Particle length-dependent titanium dioxide nanomaterials toxicity and bioactivity. Part. Fibre Toxicol. 2009, 6, 35. (17) Porter, D.; Hubbs, A.; Mercer, R.; Wu, N.; Wolfarth, M.; Sriram, K.; Leonard, S.; Battelli, L.; Schwegler-Berry, D.; Friend, S.; Andrew, M.; Chen, B.; Tsuruoka, S.; Endo, M.; Castranova, V. Mouse pulmonary dose- and time course-responses induced by exposure to multi-walled carbon nanotubes. Toxicology 2010, 269, 136−147. (18) Qian, P. Z. Sliced latin hypercube designs. J. Am. Stat. Assoc. 2012, 107, 393−399. (19) Qian, P. Z.; Wu, C. J. Sliced space-filling designs. Biometrika 2009, 96, 945−956. (20) Sacks, J.; Welch, W. J.; Mitchell, T. J.; Wynn, H. P. Design and analysis of computer experiments. Statistical Science 1989, 4, 409−423. (21) Martin, J. D.; Simpson, T. W. On the use of kriging models to approximate deterministic computer models. ASME 2004 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 2004; pp 481−492. (22) Simpson, T. W.; Peplinski, J.; Koch, P. N.; Allen, J. K.On the use of statistics in design and the implications for deterministic computer experiments.Design Theory and Methodology Conference - DTM’97, 1997; pp 14−17. (23) Kleijnen, J. P.; Van Beers, W. C. Application-driven sequential designs for simulation experiments: Kriging metamodelling. Journal of the Operational Research Society 2004, 55, 876−883.

(24) Huang, D.; Allen, T.; Notz, W.; Miller, R. Sequential kriging optimization using multiple-fidelity evaluations. Structural and Multidisciplinary Optimization 2006, 32, 369−382. (25) Ankenman, B.; Nelson, B. L.; Staum, J. Stochastic kriging for simulation metamodeling. Oper. Res. 2010, 58, 371−382. (26) Chen, X.; Zhou, Q. Sequential experimental designs for stochastic kriging. Proceedings of the 2014 Winter Simulation Conference, 2014; pp 3821−3832. (27) Porter, D. W.; et al. Differential mouse pulmonary dose and time course responses to titanium dioxide nanospheres and nanobelts. Toxicol. Sci. 2013, 131, 179−193. (28) Bonner, J. C.; et al. Interlaboratory evaluation of rodent pulmonary responses to engineered nanomaterials: the NIEHS NANO GO consortium. Environ. Health Perspect. 2013, 121, 676−682. (29) Xia, T.; et al. Interlaboratory evaluation of in vitro cytotoxicity and inflammatory responses to engineered nanomaterials: the NIEHS NANO GO consortium. Environ. Health Perspect. 2013, 121, 683−690. (30) Santner, T. J.; Williams, B. J.; Notz, W. The Design and Analysis of Computer Experiments; Springer, 2003. (31) Qian, P. Z.; Wu, H.; Wu, C. J. Gaussian process models for computer experiments with qualitative and quantitative factors. Technometrics 2008, 50, 383. (32) Sand, S.; von Rosen, D.; Eriksson, P.; Fredriksson, A.; Viberg, H.; Victorin, K.; Filipsson, A. Dose-response modeling and benchmark calculations from spontaneous behavior data on mice neonatally exposed to 2,2′,4,4′,5-pentabromodiphenyl ether. Toxicol. Sci. 2004, 81, 491−501. (33) Jones, D. R.; Schonlau, M.; Welch, W. J. Efficient global optimization of expensive black-box functions. Journal of Global Optimization 1998, 13, 455−492. (34) Martin, J. D.; Simpson, T. W. A study on the use of kriging models to approximate deterministic computer models. ASME 2003 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference, 2003; pp 567−576. (35) Hassan, R.; Cohanim, B.; De Weck, O.; Venter, G. A comparison of particle swarm optimization and the genetic algorithm. Proceedings of the 1st AIAA Multidisciplinary Design Optimization Specialist Conference, 2005; pp 18−21. (36) Law, A. M.; Kelton, W. D.; Kelton, W. D. Simulation Modeling and Analysis; McGraw-Hill: New York, 1991; Vol. 2.

3232


Kriging-Based Design of Experiments for Multi ... - ACS Publications

Recommend Documents