Population Pharmacodynamic Modeling Using the Sigmoid Emax Model: Influence of Inter-individual Variability on the Steepness of the Concentration–Effect Relationship. a Simulation Study

The relationship between the concentration of a drug and its pharmacological effect is often described by empirical mathematical models. We investigated the relationship between the steepness of the concentration–effect relationship and inter-individual variability (IIV) of the parameters of the sigmoid Emax model, using the similarity between the sigmoid Emax model and the cumulative log-normal distribution. In addition, it is investigated whether IIV in the model parameters can be estimated accurately by population modeling. Multiple data sets, consisting of 40 individuals with 4 binary observations in each individual, were simulated with varying values for the model parameters and their IIV. The data sets were analyzed using Excel Solver and NONMEM. An empirical equation (Eq. (11)) was derived describing the steepness of the population-predicted concentration–effect profile (γ*) as a function of γ and IIV in C50 and γ, and was validated for both binary and continuous data. The tested study design is not suited to estimate the IIV in C50 and γ with reasonable precision. Using a naive pooling procedure, the population estimates γ* are significantly lower than the value of γ used for simulation. The steepness of the population-predicted concentration–effect relationship (γ*) is less than that of the individuals (γ). Using γ*, the population-predicted drug effect represents the drug effect, for binary data the probability of drug effect, at a given concentration for an arbitrary individual. Supplementary Information The online version contains supplementary material available at 10.1208/s12248-020-00549-7.


INTRODUCTION
The relationship between the concentration of a drug and its pharmacological effect is often described by empirical mathematical models, as a convenient method to explore this relationship quantitatively, e.g., to predict the time course of drug effect by pharmacokinetic-pharmacodynamic modeling (1)(2)(3). This approach can be used for both continuous and binary (also denoted quantal or dichotomous) drug effects, and can be extended to combinations of drugs using response surface modeling to investigate the interaction of two or more drugs (4)(5)(6)(7).
The concentration-effect relationship is usually described by the sigmoid E max model. This model has a limited physiological and mechanistic basis, since it reflects the relationship between drug concentration and effect in the case that the drug effect is proportional to the receptor occupancy; in this case, the concentration at which the drug effect is 50% of the maximal effect (C50) equals the dissociation constant K d and the slope of the concentration-effect relationship (γ) equals 1 (see "METHODS" for details). In other cases, the sigmoid E max model should be considered as an empirical equation that describes the concentration-effect relationship sufficiently accurately, as has been shown in numerous papers over the last four decades.
An interesting characteristic of the sigmoid E max model is that the concentration-effect relationship is close to that of the cumulative log-normal distribution (8). This implies that it can be used in cases where the concentration-effect is likely to follow a cumulative log-normal distribution, e.g., in the case of binary responses, where the probability of response is modeled as a function of drug concentration. However, the similarity of the sigmoid E max model and the cumulative lognormal distribution has not been investigated in detail in pharmacological literature.
The relationship between the concentration of anesthetic drugs and binary measures of depth of anesthesia, e.g., response to a standardized stimulus, is often rather steep. However, the reported steepness varies widely between studies, even for the same stimulus and using a similar study design (9)(10)(11)(12). In some studies, the population analysis resulted in a final analysis including inter-individual variability (IIV) in C50, based on the lower value for the objective function value (10,11), but in another study no IIV in C50 could be identified (9). This results in remarkable differences in γ, which was reported as 3.46 for propofol (without IIV), 17.6 for sevoflurane + propofol (with 31-32% IIV), and 7.41 for sevoflurane (with 20% IIV). A simultaneous analysis of the data of these three studies revealed that the steepness of the concentration-effect relationship is, among other factors, dependent on the inclusion or exclusion of IIV in model parameters (12). Exclusion of IIV in C50 resulted in a lower value for steepness. This phenomenon has been investigated in two papers (3,13), demonstrating that, when data from multiple patients is naively pooled, the estimates of γ may be biased, and the 95% confidence intervals may not contain the true value. The authors stated: "We believe that estimates of γ from studies in which data from multiple patients was naively pooled must be viewed with suspicion. In this type of analysis, intra-patient variability (embodied in the parameter γ) cannot be distinguished from inter-patient variability. Accurate estimates of γ necessitate methods of analysis that take inter-patient variability into account." However, they did not provide insight in how such a study design should be chosen.
It is the aim of this paper (procedure 1) to describe quantitatively the relationship between the steepness of the concentration-effect relationship and IIV in the model parameters, and (procedure 2) to investigate whether IIV in the model parameters can be estimated by population analysis with data obtained from study designs as used in reported clinical research studies. To this purpose, several simulation studies were performed. Finally, we discuss the question whether or not IIV should be included in population modeling of binary responses.

METHODS
We describe the procedures for a continuous drug effect as well as for a binary drug effect. For convenience, we consider here the situation where only one drug is administered. Once the principle has been developed, the method can be extended for any combination of drugs, including combinations of hypnotic and opioid drugs (12).

Sigmoid E max Model
If we assume that the drug effect in an individual can be predicted from the sigmoid E max model, the drug effect P is defined: where C is the effect-site concentration, C50 is the effect-site concentration resulting in P = 0.5, and γ is a dimensionless model parameter, reflecting the steepness of the concentration-effect relationship. In the case of a binary drug effect, P reflects the probability of a drug effect at drug concentration C. Throughout this paper, it is assumed that the model parameters are log-normally distributed within the population. An example is shown in Fig. 1 (upper panel, thin lines). As a result of IIV in C50 and γ, the concentration-effect relationship is different for each individual.

Cumulative Log-Normal Distribution
If we assume that P in an individual can be predicted from the cumulative log-normal distribution function, P is defined: where Φ(z) is the cumulative normal distribution function (ranging from 0 for z = − ∞ to 1 for z = + ∞), and z is the normalized distance to the mean; for a log-normal distribution: where σ is the (dimensionless) standard deviation of the lognormal distribution, which determines the steepness of the concentration-effect relationship. An example is shown in Fig. 1 (lower panel, thin lines). As a result of IIV in C50 and σ, the concentration-effect relationship is different for each individual. Note the similarity between the profiles in both panels of Fig. 1; Fig. 2 shows that the difference in P (ΔP) is less than 0.01 (1%) over the entire scale and is minimal for P = 0.5 and at P = 0.115 and P = 0.885.

Comparison of Sigmoid E max Model and Cumulative Log-Normal Distribution
A logistic approximation of the cumulative normal distribution has been described by Bowling and colleagues (8): where z is a normally distributed random variable with mean zero and variance one, and a is a constant. The best agreement between the cumulative log-normal distribution (Eqs. (2) and (3)) and the logistic function (Eq. (4)) was found for a = 1.702 (8). For our purposes, a constant value of 1.7 was found to be sufficiently accurate (and more convenient in practice) and was used throughout this article. Equating P from Eqs. (1) and (4), it follows that Combining Eqs. (1), (3), and (5), it follows for a = 1.7 Using Eq. (6), σ (describing the steepness of the concentration-effect relationship in the cumulative lognormal distribution) can be converted to γ (describing the steepness in the sigmoid E max model), and vice versa. Note that Eq. (6) is an approximation, since both distributions are close to each other, but not identical. However, as shown in Fig. 2, the difference between the two distributions is small and likely not distinguishable when applied to clinical data.

Influence of IIV in Population Analysis
The influence of IIV in population analysis is shown in Fig.  1. Each individual (thin lines) shows a steep concentrationeffect relationship, with a different steepness due to IIV in γ (upper panel) or σ (lower panel), whereas IIV in C50 results in a shift along the concentration axis. The steepness of the concentration-effect relationship of the typical individual (dashed line) is similar to that of the individuals. However, the population-predicted concentration-effect relationship, which can be calculated by averaging P at each concentration value (solid line), is less steep than that of the individual patients. This population-predicted P represents the probability of a drug effect at a given concentration averaged across all individuals, which may be interpreted as the probability of a drug effect at a given concentration for an arbitrary individual. The steepness of the population-predicted P may be expressed by the sigmoid E max model (using symbol γ*, to discriminate from the model parameter γ) or cumulative log-normal distribution (using symbol σ*, to discriminate from the model parameter σ) and is a function of γ (or σ) and the IIV in C50 and γ (or σ).

Procedure 1: Population Predictions for σ* and γ*
Using the similarity between the probability of drug effect according to the sigmoid E max model and the cumulative log-normal distribution, an equation for the relationship between γ and IIV in C50 and γ on the population-predicted γ* was derived in a Monte Carlo simulation study.
The combined effect of γ and IIV in C50 and γ on the population-predicted γ* cannot be derived mathematically. However, the combined effect of σ and IIV in C50 and γ on σ* may be evaluated from the principle that variances sum up. According to this principle, we postulated that the variance σ* 2 may be approximated by the following expression, summing σ 2 and the variances ω C50 and ω γ (using the symbol omega in NONMEM for inter-individual variance) and the interaction of ω C50 and ω γ : where σ is defined by Eq. (6) and p i (i = 1 to 8) are constants, which were estimated using the numerical procedure described below. Note that all parameters and constants in Eq.
The population-predicted steepness γ* from Eq. (8) can be used to calculate the profile of population-predicted P.
The value of γ* calculated from Eq. (8) is lower than that of the typical individual ( Fig. 1) since it takes into account the population variability of C50 and γ.

Estimation of Constants p i in Eq. (7) by Monte Carlo Simulation
The distributions as shown in Fig. 1 cannot be used directly for population predictions. For population predictions, we need to derive the relationship between the drug concentration and P within a (large) population of individuals. This was done by estimating the constants p i in Eq. (7) with the following Monte Carlo simulation procedure: 1. A large number of individuals (here 10,000) was simulated, with randomly drawn values for C50 and γ using the sigmoid E max model. Separate sets were generated for each combination of population values for C50 (arbitrary unit with value 1 for each set), γ (values 0.5, 1, 2, 3, 5, 10, 20, 30, 50), ω C50 (0, 0.05, 0.1, 0.2, 0.3, 0.5) and ω γ (0, 0.05, 0.1, 0.2, 0.3, 0.5); the total number of generated sets was 1 × 9 × 6 × 6 = 324. 2. For each simulated individual, P was calculated over a wide range of concentrations, using Eq. (1) (sigmoid E max model) or Eq. (2) (cumulative log-normal distribution). The concentrations were equally spaced on a logarithmic scale, according to the following equation: where C i is the ith concentration (i = − 50 to 50), C50 is the typical C50 value (population mean), p is the lower range of probability levels tested (here p = 0.01, implying that P for the lowest and highest concentration is 0.01 and 0.99, respectively, and γ* is defined by Eqs. (7) and (8). 3. When simulating binary data, the actual binary status of each patient for each calculated P was simulated by drawing a random number between 0 and 1: if P was above this value, the drug effect was considered to be present (assigned a value of 1); if P was below this value, the drug effect was considered to be absent (assigned a value of 0). However, it was confirmed by simulations that this dichotomy step may be left out, since it results in the same values for P over the entire concentration range, with much better precision, avoiding loss of information in the dichotomy step. Therefore this step was skipped in all presented results. 4. At each concentration point, the population-predicted P was calculated as the sum of P for each individual (obtained by Monte Carlo simulation in step 3), divided by the number of individuals. 5. At each concentration point, the population-predicted P was also calculated using Eq. (1) with the population-predicted γ* calculated using Eqs. (7) and (8). 6. The best-fitting values of C50 and γ* were obtained by minimizing the sum of the squared differences between P calculated in steps 4 and 5 over the entire concentration range (Eq. 9), using the Solver in Excel 2019 (Microsoft, Redmond, Washington USA). 7. From the 324 combinations of γ, ω C50 and ω γ (step 1), the best fitting values of p i (i = 1 to 8) in Eq. (7) were obtained by minimizing the sum of the squared differences between the logarithms of γ* obtained in step 4 and γ* calculated from Eqs. (7) and (8), using the Solver.
To calculate the concentration values used in the simulations from Eq. (9), γ* must be known, so the constants p i in Eq. (8) must be known, whereas these simulations are intended to obtain the empirical values of p i . To solve this chicken-and-egg problem, initial values of the constants p i were obtained by a stepwise procedure: (1) p 1 and p 5 were solved from the simulation sets with varying values of γ, and fixed values ω C50 = ω γ = 0 and fixed values p 2 = p 3 = p 4 = p 6 = p 7 = p 8 = 0; (2) p 2 and p 6 were solved from the simulation sets with varying values of γ and ω C50 , and fixed value ω γ = 0 and fixed value p 3 = p 4 = p 7 = p 8 = 0; (3) p 3 and p 7 were solved from the simulation sets with varying values of γ and ω γ , and fixed value ω C50 = 0 and fixed value p 4 = p 8 = 0; (4) p 4 and p 8 were solved from the simulation sets with varying values of γ, ω C50 and ω γ ; (5) Finally, all constants p i were estimated simultaneously from all data, until the estimated values of p i were similar to that used in the preceding step.

Procedure 2: Simulation and population analysis using NONMEM
To validate the relationship between γ, ω C50 , ω γ and γ* in a situation comparable to that in reported clinical studies (9-12), a series of simulations were performed with the sigmoid E max model, followed by population analysis using NONMEM version 7.3.0 (Icon Development Solutions, Hanover, MD).
Four drug levels in each individual were chosen for each combination of model parameters in such a way, that the population-predicted values of P, calculated from Eq. (1) using the "true" population C50 and γ* from Eq. (10) were 0.10, 0.25, 0.75, and 0.90, respectively. These four values were chosen to cover a large part of the informative part of the concentration-effect profile and to guarantee that the information density was similar in each set of simulations.
At each of the four drug levels, P was calculated from the individual C50 and γ and drug level. A random number between 0 and 1 was generated; if P was above this value, the drug effect was considered to be present (assigned a value of 1); if not, the drug effect was considered absent (assigned a value of 0).
The simulated data sets were analyzed using the sigmoid E max model three times, with IIV in C50 and γ, IIV in C50 only, and without IIV, respectively. The resulting population values of C50 and γ* were evaluated by comparing the median value of 1000 replications with the "true" value of C50 and the value of γ* calculated using Eq. (11).

Calculations
The simulations and population analyses were performed using NONMEM version 7.3.0. The following code was used to control the estimation step: $ESTIMATION SIG=4 MAX=9999 METHOD=COND LAPLACE LIKELIHOOD. For the covariance matrix, the default setting was used. Other calculations were performed in Excel 2019 (Microsoft, Redmond, Washington USA).

Procedure 1: Population Predictions for σ* and γ*
To estimate the constants p i in Eq. (7), a series of simulations were performed using the Monte Carlo procedure described in the methods section. For each combination of γ, ω C50 , and ω γ , the population predictions for C50 and γ were estimated by fitting Eq. (7) to the average of 10,000 simulated probability profiles. The estimated γ* was compared to γ* obtained from Eq. (7). A selection of results is presented in Table I (complete results in supplemental table S1). From these results, the best fitting values, rounded to practical values, were p 1 = 1, p 2 = 1, p 3 = 1.25, p 4 = 2.5, p 5 = 0, p 6 = 0, p 7 = 2, p 8 = 1, resulting in the following equation: Combining Eqs. (6), (8), and (10), results in the following approximation for γ*: Note that σ* = σ and γ* = γ in the absence of IIV in C50 and γ (ω C50 = 0 and ω γ = 0). Using Eq. (11), γ* could be predicted with a mean precision of 1.2% (root mean squared error; range − 4.4 to + 5.2%) for any of the tested combinations of γ, ω C50 and ω γ . Figure 3 shows an example of the relationship between γ* and ω C50 obtained by Monte Carlo simulation (open symbols), and the model prediction using Eq. (11) (solid lines), for γ = 1, 5 and 30 and ω γ = 0. In contrast to ω C50 , the influence of ω γ on γ* is rather small; as shown in Fig. 4, γ* is hardly affected by ω γ . Interestingly, the influence of ω C50 and ω γ is rather small and comparable for γ = 1, but for γ = 30, the effect of ω C50 on γ* (and σ*) is much larger than that of ω γ . This differential effect is well predicted by Eqs. (10) and (11).
To confirm the validity of Eqs. (10) and (11), the simulated data sets were analyzed also using the cumulative log-normal distribution function (Eqs. (2) and (3)) instead of the sigmoid E max model. As shown in Table I (complete   results in supplemental table S1), the estimated values of σ were close to the values σ* calculated using Eq. (10), also with a mean precision of 1.2% (range − 4.7 to 4.5%).

Estimation of IIV in C50 and γ
To investigate whether the parameters C50 and γ as well as their IIV (ω C50 and ω γ ) could be estimated from clinical study data with acceptable precision, a series of simulations was performed using NONMEM. Each simulation consisted of 1000 runs with 40 individuals each, and 4 binary observations in each individual (see methods section for details). The results have been summarized in Table II  (complete results in supplemental table S2).
In all simulations, median C50 was very close to the true value for each tested value of γ (1, 5, 30), but 95% confidence intervals of C50 were rather wide for low values of γ. Median γ was overestimated by about 6%. For γ = 30 and IIV in C50 during simulation, γ became very high, resulting in a "near boundary" message in NONMEM. Minimization was successful in most runs, but often the covariance step was not performed, in most cases due to a "near boundary" occurrence.
When IIV was absent during simulation, the estimate of the corresponding variance became very low, resulting in a "near boundary" message in NONMEM, indicating the absence of variance in the corresponding parameter. In most cases where IIV was present during simulation, the corresponding variance could be estimated, but the estimated values were far from the variance used in the simulations. Therefore it can be concluded that the tested study design (40 individuals with 4 binary observations in each individual) is not suited to estimate IIV in C50 and γ with reasonable precision.

Estimation of IIV in C50
A similar series of simulations were performed assuming that IIV in γ was absent during simulation and estimation. The results have been summarized in Table III (complete  results in supplemental table S3). Comparing the results of  Tables II and III, it can be concluded that the influence of IIV in γ is small compared to that of IIV in C50. However, for γ = 5 and ω γ = 0.5, the estimated value for γ (18.8) is far from its true value (supplemental table S3). In addition, accuracy of estimates of ω C50 is poor. For example, for γ = 30 and ω C50 = 0.1 during simulation, the estimated ω C50 = 0.413 (Table III).

Estimation without IIV
In Table IV the results of estimation without IIV (naive pooling approach) are shown (complete results in supplemental table S4). Minimization and covariance step were successful in almost all runs, irrespective of the values of γ, ω C50 , and ω γ .
In the absence of IIV in C50, the results are broadly comparable to that in Table II, but the bias in γ is smaller (about 2% versus 6% in Table II) and confidence intervals of γ are narrower. As expected from the results in Table I, median γ decreases with increasing ω C50 and, to a much lesser degree, ω γ , and is close to γ* calculated using Eq. (11), with a mean precision of 1.9% (root mean squared error; range − 6.2 to + 2.0%) for any of the 48 tested combinations of γ, ω C50 , and ω γ .

Influence of Number of Individuals and Number of Observations per Individual
Another set of simulations was performed to investigate whether the number of individuals or the number of binary observations in each individual may affect the relationship described in Eqs. (10) and (11). For convenience, ω C50 was fixed to 0.1 and ω γ = 0 during simulation and IIV was assumed to be absent during estimation (naive pooling).

Simulations with Continuous Data
Finally, a set of simulations was performed to confirm that Eqs. (10) and (11) are valid for continuous data as well. Data sets were generated similar to that in Tables II, III, and IV, without the dichotomy step. Instead, random data error with mean zero and standard deviation (SD) 0.1 (corresponding to 10% of the full scale from 0 to 1) was added to the simulated values. IIV was assumed to be absent during estimation (naive pooling). Table I. Comparison of γ and σ estimated by fitting to the sum of 10,000 simulated probability profiles of the sigmoid E max model (γ mc ) and cumulative log-normal distribution (σ mc ) and calculated from Eqs. (11) (γ*) and (10)   The results have been summarized in supplemental table S6. The median estimates of γ were close to the values γ* calculated from Eq. (11) (Table I), with a mean precision of 2.5% (root mean squared error; range − 7.6 to + 2.3%) for any of the 48 tested combinations of γ, ω C50 and ω γ . This demonstrates that Eqs. (10) and (11) are valid for continuous data as well as for binary data.

DISCUSSION
The first aim was to describe quantitatively the relationship between the steepness of the concentration-effect relationship and IIV in the model parameters. To this purpose, we derived an empirical equation using the principle that variances sum up. According to this principle, we postulated that the variance σ* 2 may be approximated by Eq. (7), summing σ 2 and the variances ω C50 and ω γ and the interaction of ω C50 and ω γ . Using a series of Monte Carlo simulations (procedure 1), the constants p i in Eq.  (7) were estimated, resulting in Eqs. (10) and (11), describing the combined effect of γ and IIV in C50 and γ on the population estimates σ* and γ*, respectively. Note that Eq. (11) was derived using Eqs. (6) to (8), which are based on the similarity of the sigmoid E max model and the cumulative log-normal distribution. The difference between both distributions, depicted in Figs. 1 and 2, is less than 0.01 (1%) over the entire scale and likely not distinguishable when applied to clinical data. Both distributions can be interconverted using Eq. (6). Note that Eq. (6) is an approximation, since the distributions are not identical. The close resemblance of both functions implies that both functions can be chosen in population analysis, without a clear preference.
The second aim was to investigate whether IIV in the model parameters can be estimated by population modeling with data obtained from study designs as used in reported clinical research studies (9)(10)(11)(12). From our simulations (procedure 2) it can be concluded that the tested study design (40 individuals with 4 binary observations in each individual) is not suited to estimate IIV in both C50 and γ with reasonable precision. Assuming that IIV in γ is absent during estimation results in slightly more precise estimates of C50, γ, and IIV in C50. Using a naive pooling procedure, i.e., assuming absence of IIV in all parameters during estimation, results in more precise estimates. In the case that IIV is present during simulation, the estimated population estimates γ* are significantly lower than the value of γ used for simulation. The effect of model parameters and their IIV on parameter estimates has been described in a few Monte Carlo simulation studies and with binary data (3,13) and continuous data (14,15). Lu and colleagues investigated the reliability of pharmacodynamic analysis by logistic regression. Some of their findings were confirmed in our study, e.g., the good accuracy of the estimates of C50 and the minor impact of inclusion or exclusion of IIV in γ (13). They stated that, when data from multiple patients is naively pooled, the estimates of γ may be biased, and the 95% confidence intervals may not contain the true value. The authors stated in the legend to their figure 9: "A possible explanation of how the estimate of steepness of the concentration-effect relation (γ) may be biased when data from multiple patients is pooled for analysis. In this example, single data points are taken from each of nine different patients, each of whom have a steep concentration-effect relation but different values of drug concentration associated with a 50% probability of drug effect (C50) The resultant pooled concentration-effect curve appears flat (i.e., the apparent value of γ is lower than the true value." and in the conclusion they stated: "When we simulated pooled data from multiple patients (with log-normal distributions for C50 and γ), there was a larger bias in γ estimates (up to 30%), even when n was large and %CI was significantly smaller" (3). The results of Lu and colleagues are supported by the findings in the present paper. However, we believe that the qualification "bias" should not be used here. Instead of stating that the estimate of γ is biased when obtained from a pooled analysis, we propose to state that the steepness is reflected in the population estimate obtained from a pooled analysis, denoted γ* in our paper. This parameter γ* is not the steepness of the concentration-effect relationship in an individual. Instead, it describes the population-predicted value of P as a function of the concentration, and therefore it predicts the probability of drug effect in an arbitrary individual, which can be used in pharmacokinetic-pharmacodynamic information displays that allow bedside prediction of the probability of response to standardized stimulus, such as the commercially available SmartPilot® View (Draeger, Germany) and Navigator Applications Suite (GE Healthcare, USA). In contrast, estimating γ as well as IIV in C50 and γ provides a value of γ that is not directly suited for predicting drug effect in an arbitrary individual, since it does not take into account the inter-individual variability in C50 and γ. For this reason we used a naive pooling approach in the development of a triple interaction model to estimate the potency of a combination of sevoflurane, propofol, and remifentanil (12). Our simulations show that it does not seem easy to design and perform a study to estimate C50 and γ as well as their IIV. In the data of the investigated studies (9-11) we found several occasions where patients were responding to a stimulus after showing tolerance to that stimulus at a lower drug level, suggesting that γ is not very high; if γ is very high, such observation would be very unlikely. On the other hand, such observations may be a result of the preceding stimuli, and thus a methodological shortcoming of the study design. In addition, it remains to be determined whether γ is essentially identical for each patient, or that it includes significant IIV. Besides, it may be noted that the value of γ cannot be determined precisely in all cases, since its value becomes very high during the population analysis with IIV, as shown in clinical studies (9,12) as well as in the simulations in the current paper. A full analysis of optimal study design is out of the scope of the present paper. More information on sample size calculations can be found in Ogungbenro and Aarons (16) and a practical example can be found in the supplemental digital content to Weerink et al. (17,18).
In our simulations, we included 4 observations per individual. It is obvious that the location of the 4 observations per individual is important. It would have been possible to select an optimal design for each set of simulations. This would imply an extremely heavy burden on computer time and would require the implementation of an efficient procedure for optimal design. Instead, we choose a fixed design of 4 observations per individual, chosen in such a way that the "true" probability of drug effect P is 0.1, 0.25, 0.75, and 0.9 for each set of simulations.
In conclusion: An empirical equation (Eq. (11)) was derived describing the steepness of the population-predicted concentration-effect profile (γ*) as a function of γ and the IIV in C50 and γ and was validated for both binary and continuous data. The tested study design (40 individuals with 4 binary observations in each individual) is not suited to estimate the IIV in C50 and γ with reasonable precision. Using a naive pooling procedure, the population estimates γ* are significantly lower than the value of γ used for simulation. The steepness of the population-predicted concentrationeffect relationship (γ*) is less than that of the individuals (γ). Using γ*, the population-predicted drug effect represents the drug effect, for binary data the probability of a drug effect, at a given concentration for an arbitrary individual. Therefore using γ* is better suited to clinical tools, e.g., in anesthesia.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.