Mini Review

Estimating Concentration Response Function and Change-Point using Time-Course and Calibration Data

Qiang B¹*, Abdalla A², Morgan S², Hashemi P² and Peña E³

¹Department of Math and Stats, Southern Illinois University Edwardsville, USA

²Department of Chemistry and Biochemistry, USC, Columbia, SC 29208

³Department of Statistics, USC, USA

Submission: January 28, 2019; Published: March 18, 2019

*Corresponding author: B Qiang, Department of Math and Stats, Southern Illinois University Edwardsville, Edwardsville, USA

How to cite this article: Qiang B, Abdalla A, Morgan S, Hashemi P, Peña E. Estimating Concentration Response Function and Change-Point using Time- Course and Calibration Data. Biostat Biometrics Open Acc J. 2019; 9(3): 555762. DOI: 10.19080/BBOAJ.2019.09.555762

Abstract

In this paper the problem of determining the functional relationship between time and the concentration of a chemical substance is studied. An intervention drug is administered on the experimental unit from which the chemical substance (specimen) is measured. This drug is hypothesized to cause a change in the concentration level of the chemical substance a certain lag-time after the intervention. However, the concentration value could not be directly measured, but rather a surrogate response can be measured. In the time-course study, this surrogate response is measured using different electrodes which possess varied behaviors. To utilize these surrogate measurements arising from the different electrodes (sensors), a calibration study is undertaken which measures the surrogate response for the different electrodes at known concentration levels. Based on the time-course and calibration data sets, a statistical procedure to estimate the signal function and the lag-time is proposed. Simulation studies indicate that the proposed procedure is able to reasonably recover the signal function and the lag-time. The procedure is then applied to the real data sets obtained during an analytical chemistry experiment.

Keywords: Calibration study; Change-point; General linear model; Chemical concentration; Time-course experiment

Abbrevations: EU: Experimental Unit; IID: Identically Distributed; LS: Least-Squares; ML: Maximum Likelihood

Introduction

The general problem tackled in this paper is determining the time-behavior of the concentration of a chemical obtained from an experimental unit (e.u.) subjected to an intervention. Specifically, of interest is to study if a specific drug intervention causes a change in concentration. Concentration levels, however, are difficult to measure directly, especially when the e.u. could not be sacrificed, such as during a time-course study. Thus, the concentration levels are to be inferred indirectly through electric charge measurements which are stochastically related to the concentration levels. The measurement of the charge is reliant on the use of an electrode. But different electrodes possess different behaviors when measuring the charge. Thus, to be able to use different electrodes to infer the concentration levels, a calibration study using these electrodes is performed to determine the relationship between concentration and charge for each of the electrodes

The specific study that motivated the problem considered in this paper was performed in the laboratory of one of the authors (Dr. P. Hashemi) at the Department of Chemistry and Biochemistry, University of South Carolina. The chemical concentration of interest is that of serotonin, a substance that that has been implicated in affective disorders such as depression, and the e.u.’s are mice (in vivo). The intervention performed on the mice is the administration of a pharmacological agent or drug. The type of data sets obtained from the study are shown in Figures 1 & 2. Figure 1 shows charge measurements over time for five different electrodes, with panels 1 and 2 showing results for the two intervention agents or drugs: pargyline and GBR 12909, respectively. Pargyline inhibits serotonin metabolism while GBR 12909 is a dopamine reuptake inhibitor. These drugs or agents inhibit serotonin metabolism and dopamine reuptake. Figure 2 presents the calibration data for ten electrodes, together with fitted mean response curves. The five electrodes used in the time-course study for each of the two drugs came from these 10 electrodes utilized in the calibration study. For instance, the electrodes used in the time-course experiment with pargyline were the electrodes labeled 2, 3, 4, 5, and 10; whereas, for GBR 12909 the electrodes were those labeled 1, 6, 7, 8, and 9.

The goal in this study is to estimate the mean concentration response function over time based on the time course study and the calibration study. Another aspect in the time-course study is that the drug intervention was administered after an initial no-drug period. Thus, another important goal is to determine the lag-time post administering the drug after which the drug takes effect, if indeed the drug has an effect. Discipline-specific details of this motivating and focused application has been previously published [1]. We point out the statistical methodology developed in this paper and the ideas contained herein have the potential of being useful in other situations where calibrated measurements are obtained [2].

Mathematical setting

In this section we describe the postulated models for both the time-course study and the calibration study. For this purpose, let us suppose that we have a specimen from the e.u. We denote by tX the concentration level of the chemical at time t for this specimen. We denote by T the time at which the intervention is performed, that is, the time of administering the drug intervention. We shall postulate the stochastic model

where Δ is a lag-time and (⋅)g is a continuous function with If smoothness is desired, we may impose the condition that ()g⋅ is differentiable for 0.t> The time of intervention T is known, while the lag-time ,Δ together with the regression coefficients 0γ and 1,γ and the function ()g⋅ are unknown. The error is assumed to be white noise, with tε having a normal distribution with mean zero and variance 2,τ with τ also unknown. As discussed in section 1, the tX could not be directly measured. If it could be directly measured and we possess the values of ,tXs′ then we could estimate the model parameters, and also infer the change point ,T+Δ or the lag-time .Δ

To enable the determination of the chemical concentration levels ,tXs′ there are K possible electrodes that could be used for measuring the charge. At time ,t tn charge measurements using different electrodes will be taken. Thus, denoting by Y the charge and by E the electrode type, at time t tn pairs are obtained. The measured charge is affected by the chemical concentration level and the electrode type. The relationship between the charge and the concentration and electrode type is given by

where are unknown, and ∈ has a normal distribution with mean zero and variance 2,σ with σunknown. The function (⋅)I is the indicator function, taking a 1(0)-value depending on whether the argument is true(false). This is a linear model that incorporates an interaction between the concentration and the electrode type [3,4].

In order to estimate the model parameters in model (2), a calibration study is performed. See [5] for some review of calibration methods. In this study, we followed the classic inverse regression approach. First, we regress the response variable, ,Y on predictor ,X and estimate linear regression coefficients using least-squares method; thereafter, the value of an unknown ,Xis to be estimated given an observation of ,Yby subtracting the estimated intercept and dividing by the estimated slope [6-9]. Inference for the calibration parameters, is not trivial because of the presence of a normally distributed estimated slope in the denominator, which causes the inverse estimator to have infinite variance [10]. In this paper, we use Delta method to construct approximate confidence intervals for the calibration parameters. A more conservative confidence interval approach based on inverting simultaneous tolerance intervals was proposed by Scheffe [9] in literature. An alternative approach to the problem is referred to as reverse regression, when X’s are treated as the response and formally regressed on Y’s (even though the X’s are measured with negligible error). Krutchkoff [4] compared inverse and reverse regression using Monte Carlo simulations. Properties and limitations of the reverse estimators were studied by Williams [10] and Halperin [3].

In this calibration study, known levels of concentration are used, and the different electrodes are used to measure the charge. We denote by L the number of concentration levels, and these levels will be denoted by At concentration level 0,lxall K electrode types are utilized, and for each electrode type, there are m charge measurements. Thus, for the thl concentration level, there are N=LKm observations, and for the L concentration levels, there are NLKm= charge measurements. The data could be summarized as in Table 1 and pictorially depicted as in Figure 2.

With the linear models governing these charge measurements are given by

for being independent and identically distributed (IID) and having a normal distribution with mean zero and variance 2.σ These linear models could be fitted using object functions in a variety of statistical packages, such as the function lm or glm in the R statistical package [8].

Estimating parameters

Instead of using the calibration data representation presented in Table 1, for purposes of describing more concisely the estimators of parameters, we denote by vector of charge values. The design matrix is ,W which is an matrix whose ith row is

with xi the concentration level, and indicates whether the electrode type is .jWith denoting the error vector, the linear model could be written via

where the regression coefficient vector θ is

In this model, identity matrix.

The least-squares (LS) estimator of ,θ which is also the maximum likelihood (ML) estimator, is given by (see, for instance, any linear theory book such as [6])

This is an unbiased estimator of .θ The error variance 2σ is unbiasedly estimated by

where An unbiased estimator of the covariance matrix of ˆθ is

The validation of the assumptions underlying this linear model could be performed graphically or via the global method in [7].

Calibrated concentration estimators

The calibration problem will now be described. Suppose that, at a given time, we are given a charge measurement for a specific electrode: say, ()00,,ye where 0e is the electrode type and 0y is the charge on that specific electrode. More generally, suppose that M pairs of charge and electrode type are available at a given time: ()()()0101020200,,,,,,.MMyeyeye Based on these observations, what is an estimate of the concentration level, and what is an approximate ()1001%α−confidence interval for the concentration level?

Calibration based on one pair (y0,e0)

Consider first just having one pair ()00,ye of charge and electrode type observations. Denote by ()0000,xxye≡ the concentration level that led to this charge value of y0 for electrode type e 0.Based on the linear model relationship, y is a realization of the random variable

with the convention that Solving for 0,x we obtain

Of course, the error term 0∈ is not observable, but it has mean 0 and variance σ2. The expected value of given (y0,e0), is

This is the calibrated estimate of the concentration when given a charge value of 0y obtained using the electrode type 0.e By the Delta-Method [2], this will be approximately unbiased for the expectation of ()000,.xye We seek an approximation of its variance by using the delta method. For 1,2,,,jK= define the gradients

By the Delta-Method, an estimate of the variance of

An approximate ()1001%α− confidence interval for the concentration level, having observed a charge value of 0y using electrode type 0,e is given by

Wherequantile of a t-distribution with degrees-of-freedom N-2k

Calibration based on many pairs (y0,eo)

Next we consider the situation where several charge measurements are taken using possibly different electrode types. Let

Here 0my is the mth charge measurement which is obtained using electrode type 0.me From the preceding subsection, based on this particular measurement we obtain an estimate of the concentration level, given by

which has an approximate estimated variance of whose expression is obtained via (11). The M estimates of the concentration level obtained for each of the elements in ()00,ye will not be independent of each other, owing to the fact that they will all depend on ˆ.θ In fact, through the Delta-Method, we could obtain their approximate covariance matrix. However, for the sake of simplicity and practicality, we ignore the dependence among these M estimates. Under this simplified assumption, whose appropriateness will be examined later through simulation studies, we could then combine the M estimates by simply taking into account the possibly varying estimates of their variances. We recall the following well-known distribution theory result, which is easily proved using a Lagrange multiplier minimization approach

Theorem 1 Let12,,,mWWW be independent random variables with common mean μ and respective variances Among all linear combinationsso that the mean of the linear combination is still ,μ the one with the smallest variance coincides with the choice of coefficients given by

As a consequence, an approximate ()1001%α− confidence interval for the concentration level, having observed ()00,,ye is given by

Time-course study

As mentioned in section, 1 the main goal of the study is to determine the behavior of the concentration of the chemical over a period of time, where at some point T during the monitoring period, a drug intervention is performed. We suppose that over the period *0,T charge measurements using possibly different electrode types are performed at specified times The charge and electrode observations at time t1 are given by Based on these observation vectors, we find an estimate of the concentration level at time t1 to be

together with its estimate of its variance For each time ,lt we could also construct the approximate ()1001%α− confidence interval for the concentration value given by

Note, however, that these confidence intervals are not adjusted for multiplicity. On the other hand, based on the time course data and the resulting concentration estimates at each of the times of observations, given by

we could fit an appropriate regression curve that takes into account the possibly differing variability of the One of the simplest models that could be fitted to this data, under the hypothesis that upon the drug intervention at time T the concentration curve should change after a lag-time of ,Δ specifies that

and lε has mean zero and variance 2.τ Note that the 2τ in this model is not the same as the 2τ in the time-course model in (1). Here, the time of intervention T is known, whereas the lag-time Δ is not known. For a specified ,Δ this model is easily fitted using the lm command in R with the weights option enabled, which performs a weighted linear regression fit [8]. The coefficient of determination could be obtained and denoted by ()2.RΔ The coefficient of determination could be plotted with respect to possible values of .Δ The value of Δ that maximizes this coefficient of determination ()2RΔ is a plausible estimate for the lag-time .Δ Thus, a possible estimator of the lag-time is

This could be computed by fitting the model over a sequence of values of Δ and searching for the value of Δ that maximizes ()2.RΔ This computational method is the approach we implemented in this paper. The confidence interval described above at a given time point is appropriate if we only have data at one time point. However, in the time-course study, we actually obtain a series of concentration estimates for each of the time points considered. These pairs of time and concentration estimates are then used to estimate the time course model parameters. In the process of fitting the linear-parabolic curve, we could then also obtain pointwise confidence intervals for the concentration value at each of the time points, where we utilize the estimate of the error standard deviation. We surmise that the resulting point-wise confidence intervals is a better indicator of where the concentration values are. In the real-data application of the procedure in section 6, this point-wise confidence interval approach will be presented.

A real-data illustration

This section provides a more detailed description of the statistical methods performed in the data analysis in some portions of paper [1]. The data sets were obtained in the Hashemi laboratory at University of South Carolina. We limit our consideration to the intervention drug pargyline, which was hypothesized to have an impact on the concentration level of serotonin, in contrast to the intervention drug GBR 12909. Figures 1-2 present the time-course data and the fitted linear models. These figures are adapted with permission from [1] and are copyrighted from American Chemical Society. In the next section, we then present simulation studies to provide us some ideas on the properties of this procedure

A time-course study was performed leading to the data set with charge measurements over time (from 0 minutes to 121 minutes), with the intervention drug pargyline (10 mg/kg) administered at 60 minutes, via interperitoneal injection. The time-course data with pargyline as the intervention drug is depicted in the first panel of Figure 1. The different symbols correspond to the 5 different electrodes (out of 10 possible electrodes) used in measuring the charges. A calibration study was also performed where, for known concentration levels, charge measurements were obtained for each of 10 possible electrodes. The charges of each electrode are measured under three concentration values (25, 50, and 100 nM) each with 4 replicates. The data set obtained from this calibration study are the solid points in Figure 2.

A linear model, as described in (2), was fitted using the calibration data. The fitted linear models, using the calibration data, for each electrode type are depicted as the lines in Figure 2 for each of the 10 electrodes. A summary of the estimated model parameters are in Table 2. These linear models with interaction terms provide excellent fits to the observed calibration data with coefficient of determinations all above 99%. The interaction effects are all significant

Using these fitted linear models, given the charge measurements at each time point from the time course study, estimates of the concentration levels at each time point were obtained. This procedure yielded pairs of values of time and concentration estimates which are the solid points in Figure 3. A functional continuous model, as described in Section 5, was fitted to these pairs of time and concentration values using weighted regression via the lm command in the R environment. The estimate of the time-lag Δ after which the drug takes effect is ˆ3.32minutes.Δ= As mentioned earlier, this estimate was obtained by maximizing the coefficient of determination with respect to the possible values of Δ [see discussion prior to (18)]. The resulting fitted linear-parabolic model whose equation is given in (19) is shown in Figure 3.

where This plot is depicted together with the concentration estimates at each time point, which are the open circles in the plot. Point-wise confidence interval at each time point is also included in the plot. These 95% point-wise confidence intervals were constructed when the functional model was fitted to the pairs of time and concentration values using the predict. lm command in R. Details pertaining to the fitted functional relationship between time and concentration are summarized in Table 3. Observe that the coefficient for the time effect is not significantly different from which is to be expected since without drug intervention the serotonin concentration is not expected to change. We also mention that the final fitted model has an 2R equal to 95.59%, indicating an excellent 􀏐it of the linear-parabolic model for relating concentration to time for this study.

Simulation studies

In order to study the properties of the procedure described above for estimating the lag-time parameter ,Δ we performed a computer simulation study. All simulation programs were coded in R and function objects in Rsuch as lm, rnorm, etc. were utilized. The main purpose of this study was to determine if we are able to effectively estimate the lag-time parameter Δ using a time-course data and a calibration data. In the simulation, we generate the time-course data using the model

At time t and for the generated ,xt−value a charge measurement was undertaken using three electrode types: 1,2,3.m= The charge value was generated according to the model

and obtain good estimates of the model parameters using the procedure described in the preceding sections, in particular to obtain a reasonable estimate of .Δ

We performed the above-described simulation experiment using 5000 replications. The estimates of Δ are provided in the frequency histogram in Figure 4. In searching for the optimizing,Δ we used increments of 0:1, so that the possible values of ˆΔ contains only one decimal digit. The mean of these ˆs′Δis 5:03828 with a standard error of 0:4990. Thus, this mean is quite close to the true value of Δ which is 5.0. Regarding the other model parameters pertaining to the time-course portion, we summarize the results in Table 4. Histograms of these 5000 estimates for the four model parameters are depicted in Figure 5. Observe that the means of these estimates for each of the regression parameters are close to their true values

Another parameter in the time-course model is ,τ which is the standard deviation of the error component. When we fit the linear-parabolic model on the time-course and calibration data, an estimate of the error standard deviation is also obtained. However, this is not an estimate of τsince this estimates a larger value than τowing to the additional error contamination arising from the calibration study and the charge measurements at each time and for each electrode type. For instance, in the simulation study, the histogram of the 5000 estimates of the standard deviation of the error term provided in Figure 6 all exceed the true value of10.τ= The mean of these standard deviation estimates is 23.56979 with a standard deviation of 2.723542. On the basis of this modest simulation study, the statistical procedure developed in earlier sections using the calibration and time-course data sets appears successful in recovering the functional relationship between time and concentration and is also able to infer properly about the lag-time at which the drug starts to become effective

In order to further study the properties of the statistical procedure on the real data, we performed another simulation study using the estimated model parameters based on the data from the Hashemi laboratory. The time course data was generated using model (20), with 60,T= 3.3,Δ= 0.3,τ= and with ()ttZs′ IID ()0,1N variables. The values of other parameters were set to be the same as the estimated ones in Table 3. The time index took values in {}0,1,2,121. Charge measurements were obtained using five electrodes 2,3,4,5,10.m= The charge values were generated according to the model (21) with parameter values

These correspond to the estimates of the model parameters in section 6. To simulate the calibration data, we generated charge values with concentration taking values {}25,50,100x∈ for 10 electrodes, i.e. 1,2,,10,m= according to

distribution. Using the time-course data and the calibration data, it is of interest to see if we could reasonably recover the signal function in (22) and also obtain good estimates of the model parameters The resulting 5000 estimates of Δ are provided in the frequency histogram in Figure 7. The mean of these 5000 ˆs′Δ is 3:3008 with a standard error of 1:9341. The mean is indeed very close to the true value which is 3.3,Δ= indicating the possible unbiasedness of the estimator, though of course a simulation study could not establish this theoretical property. Note, however, that there were some negative estimates that arose and this could be a consequence of a wrong (we hypothesize larger values) specification of the error variances in the models. Details of the estimates for the other parameters are summarized in Table 5. Histograms of the estimates for these four model parameters are depicted in Figure 8.

Concluding Remarks

Motivated by a study in an analytical chemistry laboratory dealing with the concentration of the chemical substance serotonin, we developed in this paper a statistical procedure for estimating the functional relationship between time and concentration level of a substance, together with the change point, based on a time-course study that measures a concentration-surrogate variable (the charge in this application) and data from a calibration study. The novel aspect of this procedure is the presence of several measuring electrodes which have their unique behaviors when utilized to measure the surrogate variable (the charge in application). Simulation studies were performed to examine the properties of the proposed procedures, and the simulation results indicate that the procedure is able to reasonably recover the signal function and also the change-point in the signal function.

The procedure is also applied to the time-course data and the calibration data from the focus application obtained in the analytical chemistry experiment. Further theoretical studies and simulations will be desired to examine more carefully the procedure especially when applied to more complex settings such as when the ()g⋅ portion in the signal function in equation (1) is not parabolic or when it is non-parametrically specified. Of interest for further research are situations where the error distributions in the regression models are not normally-distributed, which may entail the use of non-parametric regression estimation methods.

Acknowledgements

E. Peña and B. Qiang acknowledge grant NIH P30GM103336-01A1 and the Biometry Core of the Center of Colon Cancer at the University of South Carolina. P. Hashemi and A. Abdalla acknowledge The University of South Carolina Start-Up Funds, The Eli Lilly Young Investigator Award in Analytical Chemistry, and NIH MH106563. We also thank Professor James Lynch and statistics graduate students Taeho Kim, Shiwen Shen, Jeff Thompson, and Lu Wang for their comments provided during regular research meetings.

References

1. Aya Abdalla, Christopher WA, Pavithra P, Srimal S, Beidi Qiang, et al. (2017) In vivo ambient serotonin measurements at carbon-fiber microelectrodes. Analytical Chemistry 89(18): 9703-9711.
2. George Casella, Roger L Berger (1990) Statistical inference. The Wadsworth & Brooks/Cole Statistics/ Probability Series. Wadsworth & Brooks/Cole Advanced Books & Software, Pacific Grove, CA, USA.
3. Max Halperin (1970) On inverse estimation in linear regression. Technometrics 12(4): 727-736.
4. Krutchkoff RG (1967) Classical and inverse regression methods of calibration. Technometrics 9(3): 425-439.
5. Irma Lavagnini, Franco Magno (2007) A statistical overview on univariate calibration, inverse regression, and detection limits: Application to gas chromatography/mass spectrometry technique. Mass Spectrometry Reviews 26: 1-18, 2007.
6. Neter J, Kutner M, Nachtsheim C, Wasserman W (1996) Applied Linear Statistical Models. Irwin.
7. Edsel A Peña, Elizabeth H Slate (2006) Global validation of linear model assumptions. J Amer Statist Assoc 101(473): 341-354.
8. R Core Team (2016) R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria.
9. Henry Scheffe (1973) A statistical theory of calibration. The Annals of Statistics 1(1): 1-37.
10. Williams EJ (1969) A note on regression methods in calibration. Technometrics 11(1): 189-192.