Jonsson EN, Wade JR and Karlsson MO Nonlinearity Detection: Advantages of Nonlinear Mixed-Effects Modeling AAPS PharmSci 2000; 2 (3) article 32 (https://www.pharmsci.org/scientificjournals/pharmsci/journal/32.html).

Nonlinearity Detection: Advantages of Nonlinear Mixed-Effects Modeling

Submitted: February 4, 2000; Accepted: August 27, 2000; Published: September 25, 2000

E. Niclas Jonsson¹, Janet R. Wade² and Mats O. Karlsson¹

¹Department of Pharmacy, Division of Biopharmaceutics and Pharmacokinetics, Uppsala University, Box 580, S-751 23 Uppsala, Sweden

²Medical Products Agency, Box 26, S-751 03, Uppsala, Sweden

Abstract

The purpose of this study was to address the question of whether the use of nonlinear mixed-effect models has an impact on the detection and characterization of nonlinear processes (pharmacokinetic and pharmacodynamic) in rich data obtained from a few subjects. Simulations were used to assess the difference between applying population analysis, ie, nonlinear mixed-effects models as implemented in NONMEM, and the standard 2-stage (STS) method as the data analysis method for detection and characterization of nonlinearities. Three situations were considered, 2 pharmacokinetic and 1 pharmacodynamic. Both the first-order (FO) and FO conditional estimation (FOCE) algorithms were used for the population analyses. Within each situation, rich data were simulated for 8 subjects at multiple dose levels. The true nonlinear model and a simpler linear model were fit to each data set using each of the STS, FO, and FOCE methods. Criteria were prespecified to determine when each data analysis method detected the true nonlinear model. For all 3 simulated situations, the application of population analysis with the FOCE algorithm enabled the detection and characterization of the true nonlinear models in at least a 4-fold lower dose level than the STS approach. For both of the pharmacokinetic settings, population analysis with the FO algorithm performed much more poorly than the STS approach. The superior detection and characterization of nonlinearities provided by population analysis with the FOCE algorithm should allow drug developers to better predict and define how a drug should be used in clinical practice in such situations.

Introduction

The presence of nonlinear (ie, dose dependent) pharmacokinetics, whether attributable to saturation in absorption, first-pass metabolism, binding, or excretion, can have significant clinical consequences¹ . Thus, the detection and characterization of any nonlinear processes, when they exist, is extremely important for the safe and rational use of a drug. The same is also true for characterization of pharmacodynamic processes. That dose-response relationships are often nonlinear is not questioned, but frequently the range of doses studied does not permit characterization of the entire dose-response relationship² , and simpler models, linear or log linear, are commonly used in practice. The problem arises if these simpler models are used to predict the response to doses beyond the range of those used to characterize the model. These predictions may result in the false belief that increasing the dose will result in a greater effect when, in actual fact, increasing the dose will result in a smaller than expected increase in beneficial effect and may also cause a disproportionate increase in observable side effects. Spigset³ showed that nonlinearity is already present at subtherapeutic doses of fluvoxamine and that this is important both in overdose situations as well as at therapeutic dose levels because a small increase in dose may result in disproportionally higher levels and increased risk of adverse drug reactions. The case where no dose-limiting toxicity prevents characterization of nonlinear pharmacokinetics must be rare, but has been reported⁴ .

Given that nonlinear processes occur, then the next step is their detection. In new drug development, the occurrence of nonlinear processes may first be evident during toxicology studies, where much higher doses than those used in the clinic are studied. This kind of knowledge may provide the drug developer with an awareness that the possibility of the presence of nonlinear processes exists, and possibly, some idea of the cause. Another route to nonlinearity detection is during Phase I/II studies, where single-dose pharmacokinetic data may not predict what is observed at steady state, or may not detect dose-disproportional changes in parameters. Once the nonlinearity is detected, then its characterization (and subsequent determination of the consequences) should be the goal.

Detection and characterization of nonlinearities may come from fitting nonlinear models or models with nonlinear components to the data. In broad terms, pharmacokinetic/pharmacodynamic models can be constructed either for data from a single individual or for data from the whole study population. The most commonly used individual approach is the so-called standard two-stage (STS) method, which involves fitting a model to each separate individual's data and then obtaining average (population) parameter estimates in a second step based on the individual results from the first step. Taking an individual approach to modeling the data can result in the use of simpler models because a single subject's data may not contain sufficient information to characterize the true model and can certainly result in different models being used for different individual subjects. This approach will, at best, make the interpretation of the results difficult, and, at worst, result in erroneous parameter estimates for the nonlinear model components⁵ . With respect to characterizing a known or suspected nonlinear system, the design of a study in which the intended data analysis method is STS usually involves more than 1 dose level. In contrast, the situation envisioned here is the study in which nonlinearity is first detected, ie, typically a single-dose study.

Population analysis through nonlinear mixed-effects modeling was originally proposed for the analysis of sparse data of the type routinely collected in a clinical setting⁶ , but it has also been shown to be useful in analyzing rich data collected in an experimental setting^7,8,9 . In a population analysis, the data from all subjects are pooled and the estimates of the population parameters are obtained in a single data analytical step. The results from the application of nonlinear mixed-effects models to the analysis of rich experimental data show that, at the very least, the same amount of information can be obtained as with the STS approach, and sometimes more. However, it is not known if nonlinear mixed-effect modeling would have an impact on the detection and characterization of nonlinear processes in rich data obtained from a few subjects who all received the same dose, ie, the kind of data that is usually present from early-phase traditional pharmacokinetic experiments. It is this issue that we sought to address. Simulations have been used to assess the difference between applying population analysis and the STS method for detection and characterization of nonlinearities in 3 settings, 2 pharmacokinetic and 1 pharmacodynamic in origin.

Materials and Methods

Data Simulation

Three different nonlinear models were simulated. The first 2 were pharmacokinetic models with different nonlinear elimination processes and the third was a pharmacodynamic dose-response model with a baseline component.

Structural Models

The first pharmacokinetic model comprised a 1-compartment model with a single saturable elimination pathway. Drug administration was a 1-time unit (see below) with long intravenous infusion (at a rate of R _o ). The drug amount in the body for individual j at time i (A _{I j} ) was simulated through integration of the following equation:

....................(1)

where C _{i j} and t _{i j} are the j ^th subject's i ^th drug concentration and time, respectively, and V _max,j and K _m,j are the j ^th subject's values for the maximum rate (amount per time) at which the drug can be eliminated and the Michaelis-Menten constant, respectively. Drug concentration-time data were computed from A _{i j} by division by the j ^th individual's volume of distribution (V _{i j} ).

The second pharmacokinetic model was again a 1-compartment model with drug administration given by a 1-time unit (see below) with long intravenous infusion. This time, however, elimination proceeded via 2 routes, one saturable and one linear. The drug amount in the body for individual j at time i was simulated through integration of the following equation:

....................(2)

where C _{i j} , t _i , V _max,j , and K _m,j are as defined previously and k _j is the j ^th subject's rate constant for the linear elimination pathway. Drug concentration-time data were computed from A _{i j} by division by the j ^th individual's volume of distribution.

The pharmacodynamic model given below was a sigmoid E_max model with a baseline response,

....................(3)

where E _{i j} and C _{i j} are the j ^th subject's i ^th true response and i ^th true drug concentration, respectively, and, Base _j , E _max,j and EC _50,j are the parameters describing j ^th subject's predose response, maximum response, and drug concentration at half maximum response, respectively. γ is the sigmoidicity factor (ie, the shape factor). The drug concentrations in the pharmacodynamic model were generated by a linear 1-compartment model where drug was administered by a 1-time unit (see below) with a long infusion.

The values employed in the simulations for each of the parameters in each of the 3 models are given in Table 1 .

Statistical Models

All individual parameter values (except γ, which was the same for all individuals) were simulated according to log-normal distributions, according to the same general equation:

....................(4)

where θ is the parameter value for the typical individual in the population, θ _j is the j ^th subject's value of the parameter, and η _θ,j is a normally distributed random value with mean zero and standard deviation ω _θ . Intersubject variability was 30% for all parameters. Residual variability was added to each of the simulated concentration/effects, as shown in the following equations for the pharmacokinetic models and pharmacodynamic model, respectively:

....................(5)

....................(6)

where C _{i j} and E _{i j} are the j ^th subject's i ^th pharmacokinetic and pharmacodynamic observation, respectively, and __{i j} and Ê_ij are the corresponding true values of these observations. The ε_{i j} are random variables with mean zero and standard deviations of 0.1 and 5, respectively.

Study Design

Data were generated for the same 8 subjects at each of 7 dose levels for both of the 2 pharmacokinetic models. The 7 dose levels were as follows: 1, 10, 35, 50, 200, 600, and 1,000 (arbitrary dose units). For each subject at each dose level, concentration time points were generated at 1, 2, 4, 6, 8, 12, 24, 36, 48, 72, 96, 120, 180, 240, 300, 600, and 1,000 arbitrary time units after the start of the infusion.

Data were also generated for the same 8 subjects at each of 7 dose levels for the pharmacodynamic model: 3, 5, 7, 10, 15, 20, and 30 (arbitrary dose units). For STS, an extra dose level of 60 was necessary to detect the nonlinearity. For each model and subject at each dose level, concentration time points were generated at 0, 0.5, 1, 2, 4, 8, 12, 24, 48, and 72 arbitrary time units after dosing.

Thirty data sets were generated for each model and dose level. The typical individual concentration time curves per each dose level for the 2 pharmacokinetic models and the typical individual concentration effect curves per dose level for the pharmacodynamic model are shown in Figure 1 .

Analysis of the Simulated Data

Models Fit to Data

For each data set (either each individual in the case of the STS analysis or all individuals for a given model and dose level in the case of the population analyses), the true model used to simulate the various data sets and a competing simpler linear model were fit to the data.

Thus, in the case of the pharmacokinetic models, the true models described in Equations 1 and 2 and a 1-compartment model with first-order elimination (Eq. 7 plus scaling by V _{i j} as described for Eq. 1 and 2) were fit to the data.

....................(7)

k _j is the j ^th subjects linear elimination constant and the other parameters are as defined previously.

For data sets generated from the pharmacodynamic model, the true model described in Eq. 3 and a power model for response, shown below, were fit to the data as follows:

....................(8)

where A _j is the j ^th subjects estimated slope and the other parameters are as defined previously.

Data Analysis

In the case of the pharmacokinetic data, simulated concentrations below 0.01 were assumed undetectable and were deleted from the simulated data sets. The resulting pharmacokinetic data sets were log-transformed prior to the analysis so that the same error structure used during the simulation could be employed in the analysis (it is not possible to use a log-normal distribution of the residual errors directly in NONMEM).

All data analysis, both the STS and the population analyses, was performed using NONMEM¹⁰ , Version 5. In the case of the population analysis, 4 algorithms were used, the first-order approximation method (FO), the FO conditional estimation method without and with centering of the conditional estimates (FOCE and FOCE CENTER, respectively) and the laplacian estimation method (LAPLACIAN). Thus, for each of the 2 pharmacokinetic and single pharmacodynamic models, results were obtained for fitting competing models by the STS method and 4 different population methods (FO, FOCE, FOCE CENTER, and LAPLACIAN).

For the STS analysis, the 2 competing models were fit to each separate individual's data. The resulting estimated parameter values were averaged and the interindividual variability was obtained by calculating the variability of the individual estimates.

For the population analysis, the 2 competing models were fit to data from all individuals for each dose level simultaneously and the typical parameter values, together with the interindividual variability associated with each parameter, were estimated directly. In addition, although γ was fixed to be the same for all individuals during simulation, the typical value was estimated during analysis.

Model Selection Criteria

For both the STS and population analyses, if both the nonlinear and linear (or power) models terminated successfully, selection between the competing models was based on the difference in the minimum objective function value (OFV) (the value of the function that NONMEM seeks to minimize). For the population analyses, the difference in OFV between hierarchical models is approximately χ² -distributed; for the STS analyses, when the amount of data per analysis object is limited, it is F-distributed. Because the models compared were not strictly hierarchical, a p-value (P < .01) was used.

For the STS to successfully identify the nonlinear model at each dose level, the following criteria had to be satisfied:

•Within each data set per dose level, 5 or more individuals had to be better described by the nonlinear model.

•More than 15 of the 30 data sets per dose level and 120 or more of the individuals per dose level (of which there are a total of 240, 8 per each of the 30 data sets) had to be better described by the nonlinear model.

For the population analysis to successfully identify the nonlinear model at each dose level, the following criterion had to be satisfied:

•More than 15 of the 30 data sets per dose level had to be better described by the nonlinear model.

For both the STS and population analyses, if one model did not terminate successfully, then the other model was judged best for that particular data set. If neither model terminated successfully, then no model was judged appropriate.

Results

Because the 3 conditional population methods (FOCE, FOCE CENTER, and LAPLACIAN) gave similar results, only the results obtained using the FOCE method will be reported and discussed.

Pharmacokinetic Model I

The number of data sets per dose level best described by the each of the 2 competing models (linear and nonlinear ) for the STS and FOCE estimation methods is shown inTable 2 . The FO method is not included in the table because it failed to identify the nonlinear model at all dose levels. The population FOCE method is able to detect the nonlinearity in the data at more than a 4-fold lower dose level than the STS method (given the dependence of the 4-fold figure on the dose levels included in the simulations).

Table 2 also displays the resulting median and range of the estimates of Km and Clint (V_max /Km). The median values for Km obtained using the population FOCE method are more accurate and are generally associated with a much narrower range of values, indicating that even when both the STS and population FOCE methods can identify the nonlinear model, the estimates obtained from the latter approach are more reliable.

Pharmacokinetic Model II

The number of data sets per dose level best described by the each of the 2 competing models (linear and nonlinear) for the 2 estimation methods is shown in Table 3 . It is again obvious that the population FOCE method is able to detect the nonlinearity in the data at a much lower dose level than the STS method (a 12-fold lower dose but, once more, given the dependence on the dose levels included in the simulations). The corresponding median and range for the resulting parameter estimates are also given in Table 3 . Again, the values presented show that the results obtained when the population FOCE method is used are more accurate and are associated with greater precision.

Pharmacodynamic Model

The number of data sets per dose level best described by the each of the 2 competing models (power and nonlinear) for the 3 estimation methods is shown in Table 4 . In this case, the STS method was unable to detect the true model over the competing power model at the original dose range studied. Thus, a higher dose level was also assessed (60 arbitrary units), but not even then was the STS method able to characterize the nonlinearity. In contrast to the pharmacokinetic models, both the FO and FOCE population methods were able to detect the true model at the same dose level. Additionally, for the FO and FOCE methods, there was a tendency for both to perform less well at the highest dose level. The reason for this may be that if the majority of the data are close to E_max , then the nonlinear process is harder to detect.

Discussion

The advantages of nonlinear mixed-effects modeling in the analysis of rich data have been expounded by Schoemaker et al¹¹ , who showed that the gain goes far beyond obtaining simple averages. There is also a growing body of evidence that the use of nonlinear mixed-effects modeling is useful when analyzing nonlinear data^12,13,14 . Using simulations, we sought to address whether application of nonlinear mixed-effects models for the analysis of rich data obtained from a few subjects can affect the detection and characterization of nonlinear processes.

We sought to classify each of the simulated studies with respect to the detection of the nonlinear model. With this aim, given the simulated replications of each study and the replication of individuals within each study, it was necessary to devise a classification criteria to judge when the estimation procedure succeeded in this task. The first level in this classification is to decide whether each unit of analysis (individual and study, for STS and nonlinear mixed-effects modeling, respectively) detected the nonlinearity. For this, we chose a p-value of .01 because the models are not strictly hierarchical. An alternative could have been Akaike's information criterion, but this leads to a less strict selection because Akaike's information criterion corresponds to a difference in OFV of 2 at 1 degree of freedom (instead of 6.6 at 1 degree of freedom, as was used herein for the population analyses). The second level of the classification is to decide on the fraction of the replications that had to identify the nonlinearity for the whole study to be classified as having detected the nonlinearity. We choose to take the liberal view that the fraction should be greater than 0.5. For the simulation replications, this is a logical choice, but for the individual level replications with the STS approach, it might be deemed too liberal. In a real situation, the required fraction of individuals supporting the nonlinear model over the linear would probably be required to be higher.

For all 3 situations studied (2 pharmacokinetic, 1 pharmacodynamic), the application of nonlinear mixed-effects models to the data analysis enabled the detection of the nonlinearity at a 4-fold lower dose level at least. The calculation of the 4-fold difference is obviously highly dependent on the dose levels selected for use in the present study. The ability of the drug developer to capitalize on the detection of nonlinearity at lower dose levels depends on when the data from the different dose levels are quantified and analyzed. In practice, it is not unusual that all data from one dose level are analyzed prior to proceeding to the clinical phase of the next (higher) dose level. The application of nonlinear mixed-effects models and detection and characterization of nonlinearities (if they are present) should enable the drug developer to proceed cautiously and perhaps to use a different dose increment than originally intended.

This paper focuses on the risk of Type II error, ie, failing to identify a true nonlinear model. The other side of the coin is to falsely detect a nonlinear model when, in fact, the data arise from a linear model -- the Type I error. In our results, this risk can be assessed by the outcome at the lowest dose levels for each of the 3 models (because the data, for all practical purposes, are linear at that point). What can be seen is that the risk for an incorrect identification of a nonlinear model is almost identical between the nonlinear mixed-effects models and STS. The only difference was with the FOCE method in the pharmacodynamic model. STS incorrectly classified one of the data sets as coming from a nonlinear model, whereas the nonlinear mixed-effects model did not.

It is of interest that the population FO method performed much more poorly than both the STS and FOCE analyses when analyzing the data from the 2 pharmacokinetic models, and it was unable to detect the nonlinearity at any of the dose levels studied. The FO method uses a first-order approximation, setting the random inter-individual variability to its expected value of 0 during computation of the population parameters. The FOCE method uses conditional estimates of the random inter-individual variability while estimating the population parameters. The dramatically poorer performance of the FO method relative to the FOCE method was not unexpected, given that the FOCE method is more appropriate in nonlinear situations and as the amount of data per subject increases¹⁵ , as is the case in the present study. Although unexpected, the poorer performance of the FO method relative to the STS method has been reported previously and was also true for the case of intensively sampled data (albeit no nonlinearity was present in the data)¹⁶ .

The problem of detection and characterization of nonlinearity in either pharmacokinetics and/or pharmacodynamics is not new and has been discussed in the literature^{14,17,18,19,20} . The 2 central points of these discussions are that the design of the study in question and the choice of method of analysis are important factors in detecting and characterizing nonlinearities. None of the previous discussions has examined the potential interaction between these 2 factors. The present results clearly demonstrate that the choice of method of analysis can affect the range of doses that need to be studied before a nonlinearity becomes detectable. The result is, perhaps, not qualitatively surprising. One of the disadvantages of the STS approach is the problem of which model to fit to the data. For example, a 1-compartment model better describes some subjects and a 2-compartment model better describes others. Nonlinear mixed-effects modeling overcomes this problem because a single model is fit to all subjects' data, and each individual subject contributes some information toward the estimation of the model parameters. The present situation is somewhat similar. The subjects with the lower clearances will exhibit the nonlinearity at lower dose levels (in our simulated cases), and nonlinear mixed-effects modeling enables the nonlinearity to be characterized across all the subjects. What is quantitatively more interesting is that the use of nonlinear mixed-effects modeling during a dose titration study could result in a more efficient and possibly less problematic trial for the investigator, because the study of doses that are inappropriately high because of the presence of nonlinearities could be avoided.

Conclusion

In conclusion, the use of nonlinear mixed-effects modeling for the analysis of data more usually analyzed by a STS approach (ie, rich data typically collected from few subjects) can provide the data analyst with a greater power to detect any nonlinearities in the data. The superior detection and characterization of nonlinearities using nonlinear mixed-effects modeling should allow drug developers and regulators to define how a drug should be used in clinical practice in these situations.

References

1. Ludden TM. Nonlinear pharmacokinetics. Clinical Implications. Clin Pharmacokinet. 1991;20:429-446. [PUBMED]

2. Dutta S, Matsumoto Y, Ebling WF. Is it possible to estimate the parameters of the sigmoid Emax model with truncated data typical of clinical studies? J Pharm Sci. 1996;85:232-239. [PUBMED]

3. Spigset O. The selective serotonin reuptake inhibitor fluvoxamine. Umeå University medical dissertations, 1997.

4. Blackledge GRP. High-dose bicalutamide monotherapy for the treatment of prostrate cancer. Urology 1996;47(suppl 1A):44-47. [PUBMED]

5. Sedman AJ, Wagner JG. Importance of the use of the appropriate pharmacokinetic model to analyse in vivo enzyme constants. J Pharmacokin Biopharm. 1974;2:161-173. [PUBMED]

6. Sheiner LB, Rosenberg B, Marathe VV. Estimation of population characteristics of pharmacokinetic parameters from routine clinical data. J Pharmacokin Biopharm. 1977;5:445-479. [PUBMED]

7. Sheiner LB, Beal SL. Evaluation of methods for estimating population pharmacokinetic parameters. ii. Biexponential model and experimental pharmacokinetic data. J Pharmacokin Biopharm. 1981;9:635-651. [PUBMED]

8. Graves DA, Chang I. Application of NONMEM to routine bioavailability data. J Pharmacokin Biopharm. 1990;18:145-160. [PUBMED]

9. Aarons O, Mandema JW, Danhof M. A population analysis of the pharmacokinetics and pharmacodynamics of midazolam in the rat. J Pharmacokin Biopharm. 1991;19:485-496. [PUBMED]

10. Beal SL, Sheiner LB. NONMEM Users Guides. NONMEM Project Group, University of California at San Francisco, 1992.

11. Schoemaker RC, Cohen AF. Estimating impossible curves using NONMEM. B J Clin Pharmacol. 1996;42:283-290. [PUBMED]

12. Hashimoto Y, Mori S, Hama N, Nakao K, Imura H, Yamaguchi M, Yasuhara M, Hori R. Nonlinear mixed effect modelling of the pharmacodynamics of natriuretic peptides in rats. J Pharmacokin Biopharm. 1993;21:281-297. [PUBMED]

13. Hashimoto Y, Odani A, Tanigawara Y, Yasuhara M, Okuno T, and Hori R. Population analysis of the dose-dependent pharmacokinetics of zonisamide in epileptic patients. Biol Pharm Bull. 1994;17:323-326. [PUBMED]

14. Hashimoto Y, Ozaki J, Koue T, Odani A, Yasuhara M, Hori R. Simulation for the analysis of distorted pharmacodynamic data. Pharm Res. 1994;11:545-548. [PUBMED]

15. Beal SL, Sheiner LB. NONMEM Users Guide - Part VII. Conditional estimation methods. University of California at San Francisco, November 1992.

16. Rodman JH, Silverstein K. Comparison the two stage (TS) and first order (FO) methods for estimation of population parameters in an intensive pharmacokinetic study. Clin Pharmacol Ther. 1990;47:151. [PUBMED]

17. Godfrey KR, Fitch WR. The deterministic identifiability of nonlinear pharmacokinetic models. J Pharmacokin Biopharm. 1984;12:177-191. [PUBMED]

18. Godfrey KR, Fitch WR. On the identification of Michaelis-Menton elimination parameters from a single dose-response curve. J Pharmacokin Biopharm. 1984;12:193-221. [PUBMED]

19. Godfrey KR, Chapman MJ, Vajda S. Identifiability and indistinguishability of nonlinear pharmacokinetic models. J Pharmacokin Biopharm. 1994;22:229-251. [PUBMED]

20. Hashimoto Y, Koue T, Otsuki Y, Yasuhara M, Hori R, Inui K. Simulation for population analysis of Michaelis-Menton kinetics. J Pharmacokin Biopharm. 1995;23:205-216. [PUBMED]