Bayesian analysis for mixed-effects model defined by stochastic differential equations

Yan Fangrong Zhang Ping Lu Tao Lin Jinguan

(1Department of Mathematics, Southeast University, Nanjing 210096, China)(2Department of Mathematics, China Pharmaceutical University, Nanjing 210009, China)(3Laboratory of Molecular Design and Drug Discovery, China Pharmaceutical University, Nanjing 210009, China)(4State Key Laboratory of Natural Medicines, China Pharmaceutical University, Nanjing 210009, China)

?

Bayesian analysis for mixed-effects model defined by stochastic differential equations

Yan Fangrong1,2,4Zhang Ping2,4Lu Tao3,4Lin Jinguan1

(1Department of Mathematics, Southeast University, Nanjing 210096, China)(2Department of Mathematics, China Pharmaceutical University, Nanjing 210009, China)(3Laboratory of Molecular Design and Drug Discovery, China Pharmaceutical University, Nanjing 210009, China)(4State Key Laboratory of Natural Medicines, China Pharmaceutical University, Nanjing 210009, China)

The nonlinear mixed-effects model with stochastic differential equations (SDEs) is used to model the population pharmacokinetic (PPK) data that are extended from ordinary differential equations (ODEs) by adding a stochastic term to the state equation. Compared with the ODEs, the SDEs can model correlated residuals which are ubiquitous in actual pharmacokinetic problems. The Bayesian estimation is provided for nonlinear mixed-effects models based on stochastic differential equations. Combining the Gibbs and the Metropolis-Hastings algorithms, the population and individual parameter values are given through the parameter posterior predictive distributions. The analysis and simulation results show that the performance of the Bayesian estimation for mixed-effects SDEs model and analysis of population pharmacokinetic data is reliable. The results suggest that the proposed method is feasible for population pharmacokinetic data.

population pharmacokinetics; mixed-effects models; stochastic differential equations; Bayesian analysis

Population pharmacokinetic (PPK) models play a pivotal role in quantitative pharmacology study. They combine the classic pharmacokinetic(PK) analysis with population statistical models, which make them able to investigate the determined factors of drug concentration in patients group quantitatively, to guide the adjustment of administration and to provide comprehensive quantitative information for a more rational and effective clinical administration regimen. In recent years, with the trend that model-based drug development (MBDD) is promoted and strongly advocated by the FDA, the population pharmacokinetic theory and practical research have been greatly developed.

Nonlinear mixed-effects modeling (NONMEM) has been proved to be a kind of widely used and effective tool for describing the pharmacokinetic (PK) and pharmacodynamic (PD) properties of drugs[1-2]. The flexibility of modeling correlation structures, with the total variation separated into inter- and intra-individual variation contributes to the rapidly growing quantity of nonlinear mixed-effects models. This modeling approach usually leads to more accurate estimation of population parameters, especially in pharmacokinetic/pharmacodynamic (PK/PD) models including random effects. However, the traditional nonlinear mixed-effects models are most frequently based on ordinary differential equations (ODEs), assuming that the inter-individual variation is independent with the measurement errors, which means that the residuals are uncorrelated, ignoring the parts of the related dynamics that we cannot predict or explain. This assumption is well-functioning to the expected distribution of assay error, but when other sources of error appear in the distribution error, a good estimator cannot be obtained in the case of ignoring the correlation which is between the structural model parameters[3]. In fact, in the actual pharmacokinetic problems, correlations between residuals are ubiquitous, the result may not be obtained exactly under the basic statistical assumption in PK modeling[4]. The NONMEM[5]is the most commonly used software for PK/PD analysis using nonlinear mixed-effects models, and it is possible to deal with correlated residual errors using an AR model[4]. The simulation results show that the introduction of a correlated residuals model may lead to better estimates of the inter-individual variations and the structural parameters.

The other approach to modeling correlated residuals is stochastic differential equations (SDEs) in the model setup. SDEs are an extension of ODEs and facilitate the ability to split the intra-individual error into two fundamentally different types: One is a serially uncorrelated measurement error, which is mainly caused by operation mistakes; the other is system error, which may be caused by mis-specified models or true random behavior of the system. This model structure includes the statistical function of the AR model, but it is more flexible with respect to specifying models for different residuals. Some PK/PD modeling books also show that the adoption of the nonlinear mixed-effects models based on SDEs rather than ODEs may be more appropriate to model intra-individual variations[6].

However, estimating parameters in nonlinear mixed-effects models with SDE is a difficult problem and not straightforward, except for simple cases. A natural approach would be likelihood inference, but the transition densities are rarely known, and thus the explicit likelihood function is usually hard to write. Actually, there is hardly any theory for SDE models at present. Kristensen et al.[7]proposed the maximum likelihood method and maximized a posteriori to estimate the parameters in PK models with SDE; however, they only focused on single subject modeling where no inter-individual variance components were estimated. Overgaard et al.[8]suggested applying the Kalman filter to approximate the likelihood function for a SDE model with a nonlinear drift term and a constant diffusion term. Torn?e et al.[9]made this algorithm to be used in NONMEM for estimating SDEs, but the NONMEM implementation cannot be used to form Kalman smoothing estimates, which is an important feature of the SDE approach, where all data is used to give optimal estimates at each sampling point. Donnet et al.[10]proposed a Bayesian inference to analyze growth curves using mixed-effects models based on stochastic differential equations and obtained good results. Struthers et al.[11]applied the Bayesian method to the multicenter AIDS cohort study. Inspired by these, we introduce the Bayesian inference to pharmacokinetic models, estimate the parameters in the nonlinear mixed-effects models based on the stochastic differential equation and help to direct the clinical test.

The present work describes the implementation of SDEs in a nonlinear mixed-effects model with parameters estimation performed by the Bayesian analysis.

1 Models and Notation

1.1 Nonlinear mixed-effects model based on ordinary differential equation

Lety=(yi)1≤i≤n=(yij)1≤i≤n,1≤j≤ni,i=1,2,…,n,j=0,1,…,nibe the true observations. Here,yijis measurement for individualiat timetij;nis the number of individuals andniis the number of measurements for individuali. Traditional nonlinear mixed-effects models usually model this process by the nonlinear mixed-effects model based on the ordinary differential equation. Formally, the classical nonlinear mixed-effects model is written as

(1)

(2)

φi～N(μ,Ω)

(3)

wherefis the possibly nonlinear function andφ=(φi)1≤i≤nare the individual parameter vectors.φiare assumed to be independently and identically normally distributed with expectationμand varianceΩ.εijare the measurement errors, and they are also assumed to be independently and identically normally distributed with null mean and varianceσ2.

1.2 Nonlinear mixed-effects model based on stochasticdifferential equation

Under the framework of ordinary differential equations, noise is only introduced through the measurement equation (see Eq.(2)). This allows the measurement noise term to absorb the whole error due to model miss-specification or true random fluctuations of the states that may ignore the correlated residuals. Considering the correlated residuals, a stochastic process is added to the state space model (see Eq.(1)). Such that the nonlinear mixed-effects model based on SDEs can be written as

dZt(φi)=F(Zt,φi,t)dt+Γ(Zt,φi,γ2)dWt

Z(t=0)=Z0(φ)

(4)

(5)

φi～N(μ,Ω)

(6)

In the nonlinear mixed-effects model based on SDEs (see Eqs.(4) to (6)), the total variance is distinguished into three fundamentally different noises: the inter-individual variabilityΩdescribing the individual difference, the system noiseγ2reflecting the random fluctuations around the corresponding dynamic model, and the measurement noiseσ2representing the uncorrelated residuals originating from measurement assay or sampling errors.

2 Bayesian Estimation

2.1 Prior specification

The Bayesian approach allows prior distributions to be incorporated with the likelihood function to evaluate the posterior distribution of the population parameters (μ,Ω),σ2,γ2in the SDE model. Thus the first step is the choice of the prior distributions. Usual diffuse prior distributions can be chosen but the resulting posterior distributions may not be proper. Therefore, we propose to use standard prior distributions suggested by Cruz-Mesia and Marshall[12]:

(7)

2.2 Posterior computation

For the ODE model (see Eqs.(1) to (3)), the Gibbs sampling algorithms are proposed in Ref.[13]. These algorithms will not be detailed here. For the SDE model, (see Eqs.(4) to (6)), we propose to use the Gibbs algorithm, including the sampling of the auxiliary random variablesφiand vectorXiof realizations of processXifor each individual at each observation time. LetX=(X1,…,Xn)∈R(n1+1)+…(nn+1)denote the vector ofnrealizations. Hence the Gibbs sampling algorithm for the SDE model is outlined as follows:

Step 1 Initialize the iteration counter of the chaink=1 and start with initial valuesσ-2(0),γ2(0),μ(0),φ(0),X(0).

Step 2 Obtainσ-2(k),γ2(k),μ(k),φ(k),X(k)fromσ-2(k-1),γ2(k-1),μ(k-1),φ(k-1),X(k-1).

Step 3 Changektok+1 and return to Step 2 until convergence is reached.

In the case when the SDE has an explicit solution, we describe an easily implemented Gibbs algorithm. When the conditional distribution of the parameter has no explicit form, we propose to approximate it using the Metropolis-Hastings random walks.

3 Model Used for Simulation

3.1 Classical one-compartment PK model

We shall simulate experiments using a one-compartment PK model, and the compound metabolism is modeled by the following single exponential attenuation modeling:

(8)

In classical mixed-effects models, the process is modeled by a deterministic function, depending on individual random parameters. Formally, the classical nonlinear mixed-effects model is written as

i=1,2,…,n;j=0,1,…,ni

(9)

For simplicity, we model on the logarithm of the noisy observed measurements.

3.2 Stochastic one-compartment PK model

We propose to introduce a stochastic term in the ODE model(see Eqs.(9)) and extend it to the SDE model with the volatility function Γ(Zt,φ,γ2)=γZt.

dZit=-kaiZitdt+γZitdWt

(10)

(11)

with solution

(12)

By discreting the SDE, the realizationXijof the SDE is Markovian:

Xi,0=logZ0

(13)

Therefore, the SDE model on the logarithm of data is defined as

(logyi0,logyi1,…,logyini)T=(Xi,0,Xi,1,…,Xi,ni)T+εi

(Xi,1,…,Xi,ni)T=

(14)

3.3 Posterior computation and inference

We propose the condition posterior distribution for the nonlinear mixed-effects model based on the SDE in details. The prior distributions of the parameters to be estimated are defined as

1/σ2～Γ(cprior,dprior)

Then the posterior distribution of the parameters is written as

(15)

where

where

(16)

(17)

(18)

(19)

Theposteriordistributionsofγ2have no explicit form, so we use the Metropolis-Hastings random walks. The Gibbs algorithm convergence is ensured by the classical convergence theorem proposed by Carlin and Louis[14], the convergence of the Metropolis-Hastings algorithm is ensured by the theorem proposed by Mengersen and Tweedie[15].

3.4 Simulations

Tab．1 Parameter estimation and 95% credibility interval

ParameterTruevalueEstimator(95%credibilityinterval)ka0.70.7132(0.6875,0.7449)σ2ka0.010.0117(0.0078,0.0183)σ20.00160.0018(0.0013,0.0020)γ20.040.0429(0.0363,0.0467)

Fig.2 Basic goodness-of-fit graph for simulation data. (a) Observed vs. population predicted values; (b) Observed vs. individual predicted value

4 Application to Metoprolol Tartrate Dissolution Data

The method is used to model the Metoprolol Tartrate dissolution data taken from Polli et al[16]. Complete data is the percentage of released drug of four types of tablet formulations of 100-mg Metoprolol Tartrate, and each experiment is repeated six times. The data were also analyzed by Let al. using the fractional differential rate equation model in 2003 and 2004[17-18]. They found that the slow dissolving test formulation was the closest to an exponential behavior, which is used here to illustrate the methods. We only use the data up to 45 min.

The percentage of Metoprolol that has not yet dissolved is modeled as Eq.(14), whereyis the measurements for experimentiat time pointj. Estimators and their standard deviation are reported in Tab.2. The estimator ofkais in agreement with comparable values found in Ref.[18]. From Tab.2, we see that there is clear difference in the performance of different methods. The SDE method performs better than the fractional differential rate equation. Fig.3 shows the percentage of Metoprolol that has not yet dissolved data and the prediction data.

Tab．2 Parameters estimations for Metoprolol Tartrate dissolution data

ParameterEstimatorbystochasticdifferentialequation(std)Estimatorbyfractionaldifferentialrateequation(std)k^a0.0232(0.0002)0.024(0.039)γ20.074(0.009)0.08(0.0054)σ2ka0.014(0.01)σ20.041(0.025)

Fig.3 Percentage of Metoprolol that has not dissolved data and the prediction data

From Fig.3, we can conclude that the predicted trajectory is in good agreement with the actual data. The credibility interval includes most of the actual data, implying the feasibility of the Bayesian approach to nonlinear mixed-effects models based on the SDE, and the Bayesian inference has a good practical application value.

5 Conclusion

We propose a Bayesian approach to nonlinear mixed models defined by stochastic differential equations. These models are the extension to the classical nonlinear mixed-effects models whose error structures are too restrictive to model some unexplained processes. The introduction of the SDE model results in a clear validation of the model which is not the case in the standard model, by adding the new stochastic term. The proposed model is proved to be useful for the pharmacokinetic modeling.An interesting area for future research is the exploration of the model with covariate. Moreover, the extension of this work to multidimensional SDEs would also be an interesting work.

[1]Aarons L, Karlsson M O, Mentré F, et al. Role of modelling and simulation in Phase I drug development[J].EuropeanJournalofPharmaceuticalSciences, 2001, 13(2): 115-122.

[2]Sheiner L B, Steimer J L. Pharmacokinetic/pharmacodynamic modeling in drug development[J].AnnualReviewsofPharmacologyandToxicology, 2000, 40:67-95.

[3]Karlsson M O, Jonsson E N, Wiltse C G, et al. Assumption testing in population pharmacokinetic models: illustrated with an analysis of Moxonidine data from congestive heart failure patients [J].JournalofPharmacokineticsandPharmacodynamics, 1998, 26(2):207-246.

[4]Karlsson M O, Beal S L, Sheiner L B. Three new residual error models for population PK/PD analyses [J].JournalofPharmacokineticsandPharmacodynamics,1995, 23(6):651-672.

[5]Kristensen N R, Madsen H, Ingwersen S H. Using stochastic differential equations for PK/PD model development [J].JournalofPharmacokineticsandPharmacodynamics, 2005, 32(1):109-141.

[6]Krishna R.Applicationsofpharmacokineticprinciplesindrugdevelopment[M]. Springer, 2004: 166-257.

[7]Kristensen N R, Madsenb H, J?rgensena S B. Parameter estimation in stochastic grey-box models[J].Automatica, 2004, 40(2):225-237.

[8]Overgaard R V, Jonsson N, Torn?e C W, et al. Nonlinear mixed-effects models with stochastic differential equations: implementation of an estimation algorithm[J].JournalofPharmacokineticsandPharmacodynamics, 2005, 32(1):85-107.

[9]Torn?e C W, Overgaard R V, Agers? H, et al. Stochastic differential equations in NONMEM: implementation, application, and comparison with ordinary differential equations[J].PharmaceuticalResearch, 2005, 22(8):1247-1258.

[10]Donnet S, Foulley J L, Samson A. Bayesian analysis of growth curves using mixed models defined by stochastic differential equations[J].Biometrics, 2010, 66(3):733-741.

[11]Struthers C A, Mcleish D L. A particular diffusion model for incomplete longitudinal data:application to the multicenter ADIS cohort study[J].Biostatistics, 2011, 12(3):493-505.

[12]De la Cruz-Mesía R, Marshall G. Nonlinear random effects models with continuous time autoregressive errors: a Bayesian approach[J].StatisticsinMedicine, 2005, 25(9):1471-1484.

[13]Carlin B P, Louis T A.BayesandempiricalBayesmethodsfordataanalysis[M]. Chapman and Hall/CRC,2000: 378-390.

[14]Carlin B P, Louis T A. Bayes and empirical Bayes methods for data analysis[J].StatisticsandComputing,1997, 7(2):153-154.

[15]Mengersen K L, Tweedie R L. Rates of convergence of the Hastings and Metropolis algorithms[J].AnnalsofStatistics, 1996, 24(1):101-121.

[16]Polli J E, Rekhi G S, Augsburger L L, et al. Methods to compare dissolution profiles and a rationale for wide dissolution specifications for metoprolol tartrate tablets[J].JournalofPharmaceuticalSciences, 1997, 86(6):690-700.

混合效應(yīng)隨機(jī)微分方程的貝葉斯分析

言方榮1,2,4張 萍2,4陸 濤3,4林金官1

(1東南大學(xué)數(shù)學(xué)系, 南京 210096) (2中國(guó)藥科大學(xué)數(shù)學(xué)系， 南京 210009) (3中國(guó)藥科大學(xué)分子設(shè)計(jì)與藥物發(fā)現(xiàn)實(shí)驗(yàn)室,南京 210009) (4中國(guó)藥科大學(xué)天然藥物活性組分與藥效國(guó)家重點(diǎn)實(shí)驗(yàn)室, 南京 210009)

采用帶有隨機(jī)微分方程的非線性混合效應(yīng)模型對(duì)群體藥物代謝動(dòng)力學(xué)數(shù)據(jù)建模,通過(guò)在狀態(tài)方程中引入隨機(jī)項(xiàng),將常微分方程擴(kuò)展到隨機(jī)微分方程．和常微分方程相比,隨機(jī)微分方程可解決群體藥物代謝動(dòng)力學(xué)模型中相關(guān)殘差問題．利用貝葉斯估計(jì)對(duì)非線性混合效應(yīng)隨機(jī)微分方程模型參數(shù)進(jìn)行估計(jì),給出群體參數(shù)及個(gè)體參數(shù)的精確后驗(yàn)分布,將Gibbs和Metropolis-Hastings算法相結(jié)合,給出參數(shù)估計(jì)值．通過(guò)計(jì)算機(jī)模擬和實(shí)例分析驗(yàn)證了方法的可靠性,結(jié)果表明利用非線性混合效應(yīng)隨機(jī)微分方程模型及貝葉斯估計(jì)方法分析群體藥物代謝動(dòng)力學(xué)數(shù)據(jù)是可行的．

群體藥物代謝動(dòng)力學(xué);混合效應(yīng)模型;隨機(jī)微分方程;貝葉斯分析

O212

s：The National Natural Science Foundation of China (No.11171065, 81130068), the Natural Science Foundation of Jiangsu Province (No.BK2011058), the Fundamental Research Funds for the Central Universities (No.JKPZ2013015).

：Yan Fangrong, Zhang Ping, Lu Tao, et al.Bayesian analysis for mixed-effects model defined by stochastic differential equations[J]．Journal of Southeast University (English Edition),2014,30(1):122-127．

10.3969/j.issn.1003-7985.2014.01.023

10.3969/j.issn.1003-7985.2014.01.023

Received 2012-10-12．

Biographies：Yan Fangrong (1978—),male, doctor; Lin Jinguan (corresponding author),male, professor, jglin@seu.edu.cn.