Treffer: The application of a Bayesian approach to the analysis of a complex, mechanistically based model.
Original Publication: New York, N.Y. : Marcel Dekker, c1991-
0 (UK279276)
Weitere Informationen
The Bayesian approach has been suggested as a suitable method in the context of mechanistic pharmacokinetic-pharmacodynamic (PK-PD) modeling, as it allows for efficient use of both data and prior knowledge regarding the drug or disease state. However, to this day, published examples of its application to real PK-PD problems have been scarce. We present an example of a fully Bayesian re-analysis of a previously published mechanistic model describing the time course of circulating neutrophils in stroke patients and healthy individuals. While priors could be established for all population parameters in the model, not all variability terms were known with any degree of precision. A sensitivity analysis around the assigned priors used was performed by testing three different sets of prior values for the population variance terms for which no data were available in the literature: "informative", "semi-informative", and "noninformative", respectively. For all variability terms, inverse gamma distributions were used. It was possible to fit the model to the data using the "informative" priors. However, when the "semi-informative" and "noninformative" priors were used, it was impossible to accomplish convergence due to severe correlations between parameters. In addition, due to the complexity of the model, the process of defining priors and running the Markov chains was very time-consuming. We conclude that the present analysis represents a first example of the fully transparent application of Bayesian methods to a complex, mechanistic PK-PD problem with real data. The approach is time-consuming, but enables us to make use of all available information from data and scientific evidence. Thereby, it shows potential both for detection of data gaps and for more reliable predictions of various outcomes and "what if" scenarios.
AN0023430080;kuo01feb.07;2019Mar13.13:14;v2.2.500
The Application of a Bayesian Approach to the Analysis of a Complex, Mechanistically Based Model.
The Bayesian approach has been suggested as a suitable method in the context of mechanistic pharmacokinetic-pharmacodynamic (PK-PD) modeling, as it allows for efficient use of both data and prior knowledge regarding the drug or disease state. However, to this day, published examples of its application to real PK-PD problems have been scarce. We present an example of a fully Bayesian re-analysis of a previously published mechanistic model describing the time course of circulating neutrophils in stroke patients and healthy individuals. While priors could be established for all population parameters in the model, not all variability terms were known with any degree of precision. A sensitivity analysis around the assigned priors used was performed by testing three different sets of prior values for the population variance terms for which no data were available in the literature: "informative", "semi-informative", and "noninformative", respectively. For all variability terms, inverse gamma distributions were used. It was possible to fit the model to the data using the "informative" priors. However, when the "semi-informative" and "noninformative" priors were used, it was impossible to accomplish convergence due to severe correlations between parameters. In addition, due to the complexity of the model, the process of defining priors and running the Markov chains was very time-consuming. We conclude that the present analysis represents a first example of the fully transparent application of Bayesian methods to a complex, mechanistic PK-PD problem with real data. The approach is time-consuming, but enables us to make use of all available information from data and scientific evidence. Thereby, it shows potential both for detection of data gaps and for more reliable predictions of various outcomes and "what if" scenarios.
Keywords: Bayesian; Hierarchical models; Mechanistic modeling; Sensitivity analysis; Summary of evidence
1. INTRODUCTION
Bayesian methods have been suggested as a valuable tool for information-gathering in the context of the mechanistic pharmacokinetic-pharmacodynamic (PK-PD) modeling normally undertaken as part of the drug development process (Ferber, [21]; Sheiner and Wakefield, [41]). This is especially the case for the design and analysis of clinical pharmacology trials. The magnitude as well as the variability of interaction between PK and PD is vital information needed for decision making during the early stages of the life cycle of a drug. In building models for these interactions, Bayesian methods make efficient use of both the data, even if such are sparse, and the prior knowledge regarding drug and/or disease state. The Bayesian approach has already been promoted as a suitable method for physiologically based pharmacokinetic (PBPK) modeling (Bois, [8]), as well as for population PK (Lunn et al., [35]). However, very few fully Bayesian studies of complex, mechanistic PK-PD models may be found in the literature. To our knowledge, the only examples of Bayesian analysis of complex, hierarchical nonlinear models in the literature concern PBPK models in the context of toxicology (Bois et al., [10][11]; Bois [6][7]; Jonsson et al. [32]; Jonsson and Johanson [29][30][31]).
In order to investigate the possibility of using the Bayesian approach in a mechanistic PK-PD setting, we present an example of a fully Bayesian re-analysis of a previously published PK-PD model. The model provided a good description of the observed neutrophil concentrations over time. However, there was not enough information present in the data to estimate all parameters of this model independently. Instead, it was necessary to fix some model parameters to reference values from the literature. In the present analysis, we use all available external information to define priors for the parameter values. We also present the summary of evidence that was derived as a basis for the choice of priors, and subsequently use the data to update these prior distributions. Thus, the main goal of this paper is to improve upon the previous model by taking all mechanistic knowledge into account in a consistent, transparent, and fully Bayesian fashion. The resulting model would be a more reliable one, and as such, more suitable for simulations and predictions. Another objective is to present the application of the method to a real PK-PD data example, thus investigating its potential strengths and weaknesses in this field for the first time.
2. METHODS
2.1. The Data
The data used in the analysis are presented in Fig. 1 and come from four studies, three healthy volunteer studies (
Graph: Figure 1 The data used in the analysis plotted against time and stratified according to assigned treatment. Each separate line represents the data from one individual.
2.2. The Previous Neutrophil Model
A mechanistically based model for the neutrophil progression in healthy volunteers and patients following acute treatment of ischemic stroke was developed previously (Jonsson et al., [28]). The reason for developing this model was that the treatment itself was believed to act by blocking the infiltration of activated neutrophils into the site of infarction. In this context, the neutrophil concentrations were believed to be useful as pharmacodynamic markers of drug effect. The model described the PK-PD in healthy and diseased individuals and had three parts. One part described the pharmacokinetics of the given candidate drug, the second part described receptor occupancy, while the third part described neutrophil progression. The present paper presents a re-analysis of the neutrophil progression part of the previous model. The neutrophil progression model (Fig. 2) was inspired by the general semimechanistic model for white-blood cell proliferation and maturation suggested by Friberg et al. ([24]), but reformulated in order to increase the mechanistic aspects of its structure. The model was described extensively in the original publication, but is also presented here in order to make the description of the summary of evidence for the definition of priors (see Appendix) easier to follow. In the manner implemented here, the model does not contain any link to the treatment, that is, no PK.
Graph: Figure 2 Graphical depiction of the structural model being analyzed.
The model consists of a stem cell compartment with auto-generation. This compartment is linked to a mitotic compartment, which, in turn, is linked to a post-mitotic compartment. Neutrophils are released from the latter compartment into a blood compartment. A reservoir compartment is linked to the post-mitotic compartment. All model parameters and symbols relating to them are presented in Table 1. The neutrophil model was mathematically described by a system of differential equations:
Graph
where
Table 1 Explanation of symbols
A number of processes were used to account for the dynamic changes in the neutrophil counts after stroke. The first process was an instantaneous increase in the neutrophil levels in plasma. This can be viewed as a change in the ratio of demarginated to marginated neutrophils (Athens et al., [3]), and was accounted for by allowing patients to have a different
Graph
where the patient-specific parameter Marg is the fractional change in the initial neutrophil levels for patients compared to healthy volunteers.
The second process was an increased release of neutrophils from the post-mitotic state to blood, which leads to a shift in the ratio of precursor neutrophils from the reservoir to the post-mitotic compartment. The increased release was assumed to decrease with time in an exponential fashion, and characterized in the model by a half-life. Therefore, a different rate constant
Graph
The third process was a feedback from the post-mitotic state to the auto-generation rate of stem cells. (Friberg et al., [24]). The feedback signal was implemented by the use of a differential equation for auto-generation in the stem cell compartment
2.3. The Adapted Neutrophil Model
Another process that accounted for the dynamic changes in the neutrophil counts after stroke was an increased residence time of neutrophil in plasma (Athens et al., [3]). This residence time would be expected to return to its pre-stroke value in due time. In the previous analysis, this return to baseline appeared to take place very slowly and patients were, therefore, allowed to have a constant, prolonged half-life throughout the observed time-period. This process was subjected to a minor modification in the present study in order to incorporate more mechanistic knowledge, as it seems reasonable to assume that the rate constant returned to the baseline value at some point in time. A parameter describing the half-life of this effect (
Graph
2.4. Statistical Model
The statistical model has two major components: the individual level and the population level. For each of the
At the population level, the inter-individual variability is described by assuming that each vector of
2.5. Definition of Priors for the Full Model
In the previous modeling of the same data (Jonsson et al., [28]), it was found that not all parameters were identifiable from a frequentist perspective. The mean transit times through the mitotic (MRT<subs>2</subs>) and post-mitotic (MRT<subs>3</subs>) compartments were all fixed to 120 and 168 h, respectively (Walker and Willemze, [47]). The initial ratio between the levels in the post-mitotic and the reservoir (
In cases when we want to incorporate a considerable amount of information, such as here, the analyst has several choices (Spiegelhalter et al., [45]). In our opinion, the most attractive ones in the context of mechanistic PK-PD is either eliciting the opinion of experts or making a summary of available evidence. There is a considerable amount of literature on the subject of prior elicitation. However, as we, as authors, feel comfortable defining the priors ourselves, we refer the interested reader to a review (Chaloner, [13]), rather than discussing the issue of prior elicitation further.
Thus, priors for all model parameters were defined with regard to the external sources of information available in the scientific literature. Multivariate lognormal distributions were chosen for μ and σ<sups>2</sups>, while multivariate inverse gamma distributions were chosen for Σ. While defining the priors, we tried to approach the reference literature in an as objective fashion as possible, while always giving priority to values derived from in vivo studies in humans. When cross-species or in vitro extrapolations were necessary, they were made, but compensated for by an increased uncertainty around the parameters in question. A full description of the summary of evidence, including assumptions and considerations made, can be found in the Appendix. The prior distributions are listed in Table 2.
Table 2 Prior estimates of the parameters of the neutrophil model
Ideally, from a sensitivity analysis point of view (see below), we would want to also do testing around these priors that are established from the literature, as we have no idea whether they are biased or not. However, given model complexity and run times, we limited the scope of the present analysis to be performed while disregarding possible bias. Granted, a study of the effect of making a different summary of the available evidence, thus producing a different prior for elements of μ would be an interesting aspect of the analysis and model system to be investigated in future studies.
Another issue is that of correlation between parameters, both on the population and individual level. As multivariate priors allow for prior correlation among parameters, it would theoretically be possible to incorporate correlation in the priors to some degree. The only prior information available on between-parameter correlation is that provided by the covariance matrices that were estimated in the previous NONMEM runs. However, these runs were performed with a statistical model that did not allow for variability in all model parameters. Meanwhile, there is bound to be considerable information regarding parameter correlations in our data. Therefore, we do not account for any explicit prior correlations, preferring to let the correlations manifest themselves in the multivariate posterior. Parameter correlations are then easily monitored by calculating a correlation matrix for the posterior distribution.
In addition, the observations in individual patients are associated with some degree of intrasubject correlation. This correlation is no problem, as we use individual parameter estimates in conjunction with the differential equations in order to derive individually predicted concentration-time profiles.
2.6. Bayesian Model Fitting
The complicated high-dimensional joint posterior distribution was summarized by random draws using Metropolis-Hastings Monte Carlo sampling (Smith, [42]). Five independent Markov chains were run to approximate equilibrium. Posterior distributions were obtained by running the MCMC simulations further. When we subsequently derive posterior predictions, we do not assume that these parameters are independent, but sample directly from the posterior distributions.
The MCSim software (Bois and Maszle, [9]), version 5 beta, was used throughout the study. When MCSim is used, the structural model is defined, using a C-like grammar, in a model description file. The prior distributions and the statistical model are defined in an input file, which includes the information present in the clinical data. These files are then used in conjunction to update the priors using Markov chain Monte Carlo simulation. MCSim is available for free download at http://toxi.ineris.fr/activities/toxicologie%5fquantitative/mcsim/mcsim.php.
2.7. Sensitivity Analysis of Assigned Priors
In a standard, that is, frequentist, setting, a sensitivity analysis of a mechanistic PK-PD model would encompass an assessment of the extent to which each parameter affects the predictions according to some likelihood-based criterion. The least influential parameters would be fixed to their reference values, and it would then be deemed acceptable to estimate the remaining model parameters using the usual likelihood-based methods. However, when this approach is used, the fact that not all reference values are known with an equal amount of precision is largely neglected. This may not matter much for the results in cases where many model parameters are relatively well-known, and only relatively few may be regarded as very uncertain. However, when there is non-negligible uncertainty around several model parameters, the sensitivity analysis may lead to potentially misleading conclusions (Spear and Bois, [44]). In contrast, the use of a fully Bayesian approach already allows for incorporation of different degrees of uncertainty among model parameters, regardless of the magnitude of uncertainty and the distribution of this uncertainty among model parameters. This undermines the fruitfulness of a sensitivity analysis, in the frequentist interpretation of the term, that is, investigating the impact of changing one parameter while conditioning on all the others as being fixed. In one way, the MCMC process in itself may be seen as a more sophisticated and thorough investigation of the impact of changing one parameter at a time, as for every step of the Markov chain, values of a certain parameter will be drawn randomly and the likelihood calculated while conditioning on the current values of the other model parameters.
However, the results of a Bayesian analysis will always be dependent on the prior in some way or another and the definition of priors is always a subjective procedure. There is always a need for some kind of check about to what extent the conclusions are sensitive to the prior. In a Bayesian setting, the term "sensitivity analysis" is commonly interpreted as alluding to a comparison of the posterior distributions derived using informative and noninformative (so-called "reference") priors. Such a sensitivity analysis is rather straightforward for simple models (Spiegelhalter et al., [45]), but has rarely been performed in any formalized manner for complex models, such as the one in the present example.
Our adapted approach to the sensitivity analysis of the present Bayesian model is to test three different priors. In all priors, the same literature-based values were used for the population means and the variability terms for which literature data were available. The unknown variability terms were treated differently in the three priors:
Another issue that is important in the context of sensitivity analysis is the prior distribution of the error term σ<sups>2</sups>. It would seem reasonable to use a flexible prior for this parameter, as one could expect all the information on the difference between predictions and observation to be readily available during the model optimization procedure. Inverse gamma distributions have been suggested for the prior distribution of σ<sups>2</sups> in this kind of setting (Lunn and Aarons, [34]; Lunn et al. [35]). With this in mind, an inverse gamma distribution was defined, with a mean corresponding to an error of 5% on the log scale, and a relatively large standard deviation of 25%. However, when this uncertainty distribution of σ<sups>2</sups> was used in conjunction with the "informative" prior, the Markov chains did not converge, but rather showed a tendency to become more disperse as the chains progressed. This phenomenon may seem bewildering at first, but upon closer thought, one may realize that an inverse gamma prior is flexible enough to accommodate even very large residual errors. If, as in this case, the model is complex and the uncertainties around some parameters are considerable, this outcome is perfectly understandable, if not desirable. We then decided that large values of σ<sups>2</sups> are not really reflective of the degree of precision that is to be expected, given the earlier, frequentist-based modeling of the same data, and an even more flexible space for the model parameters than previously. Rather, we would expect σ<sups>2</sups> to be at most of a similar magnitude and shape as in the previous modeling of the same data (Jonsson et al., [28]). Since we did not believe much in large values of σ<sups>2</sups>, we decided to exchange the undesired flexibility of an inverse gamma prior for a more robust and narrow log-normal distribution. Thus, the error term that was estimated by Jonsson and coworkers was adapted directly as a lognormal prior of 0.0475 with a standard deviation of 0.0065. When a log-normal σ<sups>2</sups> prior based on the results from the previous model was used instead, convergence could be achieved.
2.8. Convergence Analysis
The common convergence criterion for multiple chains (Gelman and Rubin, [27]), as implemented in BOA (Smith, [43]), was used as convergence criterion, as well as visual inspection of the chains.
For every set of runs, five Markov chains were run. When a set of Markov chains had reached approximate convergence, each chain was then run for 25,000 more iterations. Every 25th iteration was printed to the output file, and thus a total of 1000 posterior parameter vectors were produced from each of the five chains.
3. RESULTS
Using the "informative" prior, five MCMC chains, using differing random seeds as starting points, were run for 150,000 iterations. One of the chains required a restart with a different seed due to slow mixing. Each chain was then run for a further 50,000 iterations while writing to file. By analyzing the first 25,000 written iterations from each chain, it was concluded that approximate convergence had been achieved. In total, the 200,000 iterations took approximately 70 h each to run on a medium-sized cluster running RedHat Linux version 9.0 with OpenMosix patch(2.4.22-3) on nodes with dual Intel Xeon processors between 2.0 GHz and 3.2 GHz and some older nodes with Intel Pentium III (Coppermine) between 0.9 and 1.4 GHz. The nodes had between 512 and 1024 Mb of Random Access Memory. The cards were interconnected with two Bonded 1 Gbit network cards. The written output from the last 25,000 written iterations was pooled and the resulting posterior distributions are presented in Table 3. From Table 3, it may be gathered that the posterior variability terms were generally a magnitude or two smaller than in the prior. In Table 4, the posterior correlation matrix is presented. The posterior distributions for most parameters show no sign of correlation. However, some of the population parameters were correlated.
Graph: Figure 3 Observed data versus means of predictions derived using posterior simulations and disregarding interindividual variability. The line of unity is also indicated.
Graph: Figure 4 Observed data versus means of predictions derived using posterior simulations where individual parameter estimates are used. The line of unity is also indicated.
Graph: Figure 5 Observed neutrophil levels and predictive 95% intervals versus time in a selection of patients.
Graph: Figure 6 Observed neutrophil levels and predictive 95% intervals versus time in a selection of healthy individuals.
Table 3 Posterior estimates of the population parameters of the neutrophil model
Table 4 Correlation matrix for the posterior distribution of the population parameters derived using an informative prior
When the "semi-informative" and "noninformative" priors were used, several correlations between parameters, most importantly
4. DISCUSSION AND PERSPECTIVES
In the present paper, we used all available external information (see Appendix) to define priors for the physiological parameters in a previously published complex, mechanistic PD model. We also presented the summary of evidence that was derived as a basis for the choice of priors, and used the previous data to update these prior distributions. All mechanistic knowledge was thus taken into account in a consistent and transparent fashion. The present study is a first attempt at applying a fully Bayesian methodology to a real problem concerning a hierarchical mechanistic PK-PD model. The results presented in the Table 3, while relying on certain assumptions regarding variability, add considerable mechanistic weight to the value of the results of the analysis, compared to the previous analysis of the same data. The methodology has potential advantages within the field of drug development. The Bayesian model would provide a more reliable basis for subsequent simulations and predictions, but such exercises lie beyond the scope of the present study, and would require added attention to the PK of the treatment, which was not included in this part of the original model (see Methods).
While we found that even the use of an informative prior at the population level in conjunction with a lognormal uncertainty distribution for the residual error was not sufficient to stabilize the joint posterior distribution entirely, in a "real" drug development setting these results could as well be interpreted as important information on data gaps. We envision that this type of analysis would typically be done in the early phases of development, and the results may then influence the design of subsequent studies. The project time spent on consulting the reference literature in order to derive prior parameter distributions, and running the multiple MCMC runs may be time well spent, if the payoff is a truly mechanistic perspective and an optimal use of information when planning future studies. Furthermore, the extra time needed for this type of analysis is likely to decrease considerably in due course. Computer power increases daily so the run time limitations we have experienced should soon disappear. The knowledge on a certain disease system or therapy area needed to build these models are already there among the members of the project team, However, focus is needed in turning the general qualitative understanding of any system through a process of critical assessment into one that is truly quantitative and useable as a prior. We think that the greater focus on modeling and simulation within drug development process will enable a transition to more Bayesian thinking and that teams will learn a lot about how well they know the disease system of interest on the way to defining the priors. In cases where there is enough information in the data for a stable joint posterior parameter distribution, the advantages of working in a truly mechanistic setting, and making all subjective assumptions transparent to the reader of the report or publication, are obvious, and the present study may be seen as merely another step along that way.
We did not achieve the goal of propagating data, prior information, and model into a fully Bayesian final model using the "semi-" or "noninformative" priors. The reason why the Bayesian approach failed in this respect may be due to methodological reasons. We adapted the generally accepted, and commonly used, choice of inverse gamma distributions for the uncertainty around the vector of population variances ∑. The inverse gamma distributions are recommended in this setting because of their flexible properties. However, it has recently been suggested (Gelman, [26]) that in cases where there is little evidence of variability in the data, inverse gamma distributions tend to not exhibit the desired flexible properties, leading to a posterior that is in fact dominated by the prior distributions. Suggested solutions to this problem include replacing the inverse gamma distributions by, for example, half-Cauchy (Gelman, [26]) or half-Normal (Spiegelhalter et al., [45]) distributions. This would definitely be an interesting direction for future work. Another possible approach would be to consider using so-called "frequentist priors" to account for differing values of prior uncertainty (Olsson Gisleskog et al., [39]), instead of adapting a Bayesian approach proper. However, this approach is relatively new, and requires the use of undocumented features in the software NONMEM. This method does not make full use of the advantages inherent to a Bayesian analysis in terms of the treatment of parameter covariances or correlations. Furthermore, in the manner that this method is presently implemented in NONMEM, it is not very flexible in terms of options for the shape of the prior uncertainties, which, as discussed in the previous section, may be of major importance in this setting.
The convergence problems we encountered (see Results) may be a signal indicating that the original structural model is an oversimplified representation of the complex mechanistic process of neutrophil maturation and does not account properly for the processes that take place in the human subjects that gave rise to these particular data, that is, model misspecification. Another possible conclusion could be that there is insufficient information in either prior or data. This conclusion could be interpreted as a signal to collect more data on the model parameters in question. This would very well have been a useful piece of information if the modeling work had been performed as part of a larger development scheme, as it efficiently points us to information gaps in the database. If such data gaps are not properly identified, they may undermine the reliability of subsequent simulations and predictions.
We conclude that the present analysis represents a first example of the fully transparent application of Bayesian methods to a complex, hierarchical mechanistic PK-PD problem with real data. This approach may very well prove useful in due time, but further work still needs to be performed in order to develop a suitable methodology for making use of its advantages. Presently, the approach is time-consuming, but enables us to make efficient use of all available information from data and scientific evidence. Thereby, it shows potential both for detection of data gaps and for more reliable predictions of various outcomes and "what if"-scenarios. This adds weight to subsequent simulation of future studies. The present study may be seen as a first step along the way to a standardized strategy for using these methods in drug development.
5. APPENDIX: SUMMARY OF EVIDENCE FOR ASSIGNMENT OF PRIORS TO MODEL PARAMETERS
Prior distributions were assigned to model parameters while taking the reference literature into consideration. Priority was always given to parameter estimates that could be based on data from humans and on in vivo rather than in vitro studies, but in some cases, animal and/or in vitro data were used in lack of other sources of information. As a rule of thumb, the uncertainty was doubled when animal data was used, and "large" uncertainties were assigned in cases where there is reason to believe that the reference estimate is not all that reliable or that the link to the estimand of interest is weak.
For some model parameters, namely variance parameters ω<subs>
5.1. k max of Exponential Function for Change in MRT 3 in Diseased Individuals (k max)
The duration of the post-mitotic phase has been estimated in healthy and diseased individuals (Fliedner et al., [23]). One of the diseased individuals in the study was observed in a healthy state, during a slight urinary tract infection and during pneumonia. The transit time for this subject was reduced by 24 and 48 h during the respective disease states. In the present model, the reduction in transit time was defined as the maximum reduction of transit time in absolute terms, with a corresponding half-life (
5.2. Stem Cell Production Rate (kn 1)
The stem cell production rate "has been poorly characterized from a kinetic point of view" (Cronkite and Fliedner, [14]), and this number has been represented by question marks in tables of suggested rate constants (Cronkite and Fliedner, [14]; Cronkite and Vincent [15]).
Thus, instead of quantifying the rate of stem cell production directly, we adapt an approach based on a semimechanistic pharmacodynamic model for neutrophil progression in conjunction with chemotherapy-induced myelosuppression (Friberg et al., [24]). In the model presented by Friberg and coworkers, the same kinetic constant was used for egresses from all compartments. The kinetic constant was estimated from the data and compared to an earlier half time value of 6.7 (Cartwright et al., [12]). This is the exact estimate that is used as a basis for MRT in plasma (MRT<subs>4</subs>) in the present model. However, we opted for a prior for MRT<subs>4</subs> of 9.85 h, with a log-normal uncertainty of 1.227. If we take the inverse of MRT<subs>4</subs>, we are left with a prior mean of 0.102. Since the original model (Jonsson et al., [28]) contained a hard-coded correlation with MRT<subs>4</subs>, the prior mean of this parameter was revised to the prior mean of MRT<subs>4</subs> divided by 0.102, producing a revised prior of 96.6.
As this value has a rather tenuous link to the actual quantity of interest, we assign a large CV of 100%.
5.3. Variability in Stem Cell Production Rate (ω kn1)
The structural model being used in the present analysis (Jonsson et al., [28]) hinges on the fact that the stem cell production rate correlates to the elimination rate from blood in healthy subjects. In the previous model, this was accomplished by using the same variability term for these two parameters.
For the variability in stem cell production rate, we have virtually no information (see section on mean stem cell production rate). However, for the variability in the elimination rate, we have no less than eight studies (Alexanian and Donohue, [1]; Cartwright et al. [12]; Dancey et al. [17]; Dresch et al. [20]; Galbraith et al. [25]; McMillan and Scott [38]; Mauer et al. [37]; Uchida [46]) from which variability may be estimated
Five of these estimates lie in the 0.22–0.27 range (quantified as relative standard deviations), while two estimates are outliers at 0.389 (Dresch et al., [20]), 0.311 (McMillan and Scott, [38]) and 0.115 (Galbraith et al., [25]). A weighted average of these eight estimates produces a population CV of 0.233, with a SD of 0.058 (CV 0.249). For the population variances, the estimate is 0.0544, with a variance of 0.00337. The alpha term for an inverse gamma distribution is then 2 + 0.0544<sups>2</sups>/0.00337 = 2.88, while the beta term is 0.0544*(1 + 0.0544<sups>2</sups>/0.00337) = 0.933.
5.4. Ratio of Number of Neutrophils in Post-Mitotic to Reservoir Compartment (k p)
In an extensive review (Osgood and Seaman, [40]), the weighted means of 110 cases from the literature were taken into account and reference values were suggested. The fraction of band form cells among the post-mitotic pool was 0.527, with a CV of 0.191.
In a later study (
Using the raw data from another study (Donohue et al., [19]), the mean of the estimates of the fraction of band form cells in the post-mitotic pool can be calculated to be 0.428, with a CV of 25.2%.
Judging from data from another study (Lala et al., [33]), the mean of the estimates of the proportion of band form cells in the post-mitotic pool was 0.348. The authors reported a "standard error" for the distribution of cell types. If this standard error is multiplied by the root of n, we may derive an estimated population CV of 22.1%.
A weighted mean of these four estimates gives 0.497, with a CV of 0.0916.
If 49.7% of the cells are in the reservoir compartment, the rate constant for the influx to the reserve compartment is then 0.497/(1 − 0.497) = 0.988. Efflux from the reserve compartment is fixed to 1, as the parameter is defined as a ratio.
5.5. Variability in Ratio of Number of Neutrophils in Post-Mitotic to Reservoir Compartment (...
A weighted mean of the four estimated population variabilities listed in previous section gives a population CV of 0.210, with a standard deviation of 0.0252. The alpha term for an inverse gamma distribution is then 5.063, while the beta term is 0.179.
5.6. Margination Effect (Marg)
De-margination of cells has been studied after epinephrine injection (
Since it is not at all certain that the de-margination observed after epinephrine injection is directly comparable in magnitude to the one that occurs after stroke, we set a high prior uncertainty of 100%.
5.7. MRT in Mitotic Phase (MRT 2)
According to a review article (Mary, [36]), there are a number of studies from which this transition time may be estimated. Unfortunately, few of these references include estimates of the residence time in the myeloblast stage, and, more importantly, all concern in vitro studies only. Therefore, we choose to disregard these studies, and instead adapt the view given in another review (Walker and Willemze, [47]), which emphasizes that a composite MRT of 11.4 (273.6 h) days, with a CI of 9.2–13.6 (
We propose a prior for the post-mitotic transit time of 135.6. If we subtract 135.6 from this estimated composite value of 264.8, we get 129.2.
Thus, we assign a lognormal prior of 129.2, with a CV set at 0.0692.
5.8. Variability in MRT in Mitotic Phase (ω MRT2)
We suggest the use of the reported CI of 9.2–13.6 (Cartwright et al., [12]) as basis for the prior estimate of this population variability term. We believe that the CI reported by Cartwright and coworkers really reflects variability between study subjects, rather than uncertainty, as there is no other source of uncertainty around the mean than interindividual variability, given the design of the study. If we assume that the confidence interval was defined by Cartwright and coworkers using the normal approximation as 1.96 times the standard deviation, the SD would be 1.122, and the population CV is then 0.098. Another estimate of population CV is 0.168 (Uchida, [46]). A weighted mean of these two estimates produces the value of 0.111, with a variance of 0.000792. The alpha and beta term for an inverse gamma distribution are then 2.195 and 0.0148, respectively.
5.9. MRT in Post-Mitotic Phase (MRT 3)
A transit time of 6.6 days (
According to a review article (Mary, [36]), there are several other studies from which additional estimates of this transition time may be derived. Since these studies are case studies with 1–3 individuals, concern mainly diseased subjects, or are only based on in vitro studies, we choose to disregard these results in favor of these two oft-cited studies. Thus, we calculate a weighted mean of the two estimates (Dancey et al., [17]; Fliedner et al. [23]), and get 140.9 (CV 16.4%).
5.10. Variability in MRT in Post-Mitotic Phase (ω MRT3)
A weighted mean of the two population CVs (Dancey et al., [17]; Fliedner et al. [23]) produces the estimate 0.0845 (SD 0.0770, variance 0.00592). The alpha and beta term for an inverse gamma distribution are then 2.00859 and 0.00720, respectively.
5.11. MRT of Neutrophils in Plasma (MRT 4)
The half-life for neutrophils in plasma has been reported as 6.6 with a standard deviation of 1.4 (
Other half-lives have been found using the <sups>32</sups>P-DFP label: 3.8 h (CV 23.7%,
In a later publication (Dancey et al., [17]), the same half time was studied using two different labels, <sups>32</sups>P-DFP and <sups>3</sups>Htdr (
<sups>32</sups>P-DFP has been found to be a less reliable label than <sups>3</sups>Htdr (Deubelbeiss et al., [18]). However, the estimate derived using the <sups>3</sups>Htdr label (Dancey et al., [17]) is very close to the values reported by the most extensive <sups>32</sups>P-DFP studies (Cartwright et al., [12]; Galbraith et al. [25]; Mauer et al. [37]). In addition, the interstudy variability among all the <sups>32</sups>P-DFP studies is larger than the interlabel variability in the single comparative study. With this in mind, we choose to accept the literature estimates at face value, while giving <sups>3</sups>Htdr-derived estimates priority over the ones derived using <sups>32</sups>P-DFP, where applicable. The weighted mean of these eight estimates is 6.83 and the CV is 20.4%.
It must be kept in mind that these values are half-life values, and need to be divided by the log of 2 to derive corresponding MRT values. This transformation produces the value of 9.85 h.
In a review article (Finch et al., [22]), a MRT value of 0.3 days (7.2 h) was tabulated, while citing the references above. A MRT of 0.3 would correspond to a half time of approximately 5 h. However, it is likely that MRT and half time were confused in this publication and this estimate is disregarded.
There are studies in the literature where a third label, <sups>51</sups>Cr, has been used, either alone or in comparison with the <sups>32</sups>P-DFP label. However, this label has been criticized for being very unpredictable and in conflict with the results derived using other labels (Finch et al., [22]). Thus, the <sups>51</sups>Cr-derived half-life estimates in the literature were disregarded during the elicitation process.
5.12. Variability in MRT of Neutrophils in Plasma (ω MRT4)
The structural model (Jonsson et al., [28]) requires that the stem cell production rate correlates to the elimination rate from blood in healthy subjects. This is accomplished by using the same variability term for these two parameters, and the reader is thus referred to the section on variability in stem cell production rate (ω<subs>
5.13. Change in MRT 4 in Patients (Mstep)
There is very little prior information for this parameter. Because of lack of other information, we have resorted to studies of the change in neutrophil plasma half-life estimated in conjunction with other diseases states than stroke.
A relative increase in MRT of 0.517 has been reported (Athens et al., [2]), but no standard deviation, only a range, in five prednisone-treated patients.
A relative increase of 0.463 (SD 0.486) (
Neutrophil kinetics in conjunction with inflammation have also been studied (Bishop et al., [5]). When the neutrophil-time curves were analyzed by curve-fitting, the egress rate constants showed a mean decrease of 74% (
If we compute the weighted mean, based on the number of subjects in the two studies where half-lives were measured directly, we get 0.385 and a CV of 0.252. 0ur proposed prior is thus lognormal, with a mean of 0.383.
There are also some studies where the mean plasma half time among small populations of infected patients have been studied (Galbraith et al., [25]; McMillan and Scott [38]). These results have been disregarded, as no pre-inflammatory half times were reported in these studies, and variability in half times generally tends to be fairly large, even in the absence of infection.
5.14. Variability in Change in MRT 4 in Patients (ω Mstep)
The sources of information are the same as for the mean change in MRT<subs>4</subs>.
Among patients with various infections (
If we choose to regard 105 and 129% as the best available estimates of the mean increase in diseased subjects, and compute the weighted mean and standard deviation of these estimates, a population standard deviation of 118%, with a SD of 0.121 and a variance of 0.0147, is produced. The alpha and beta terms for an inverse gamma distribution are then 133.745 and 184.749, respectively.
5.15. Half time of Exponential Function for Change in MRT 3 in Diseased Individuals (t half)
The apparent entry rate of granulocytes into the blood after sterile inflammation has been studied in five dogs (Cronkite et al., [16]). The apparent entry rates were reported only in aggregate form. If the average of the entry rates is plotted longitudinally, it may be noted that the apparent entry rate increases very quickly to a peak value during the first 4–6 h after inflammation. The decrease follows a biphasic pattern with approximate alpha and beta half-lives of roughly 5 and 35 h, respectively. The beta phase begins at roughly 14 h, and thus, given the time frame of data collection in the studies we want to describe, the beta half-life seems to be the relevant one for this particular parameter. We assign a prior of 35, with a lognormal uncertainty of 100%, as this value is only based on a single study in dogs.
5.16. Half Time of Exponential Mstep Function (t half Mstep)
We have no useful information on the return of the elimination half-life to pre-stroke values. It seems reasonable to assume that the effect of the disease on transit time should subside with time. In order to derive a proper estimate the rate of recovery from the accelerated transition state, several consecutive transition studies would have to be performed in the same individual before, during and after a certain disease state. That such studies have not been undertaken may not be surprising, but unfortunate for us. Still, differences in disappearance time of labeled cells in the blood stream between healthy and diseased patients are easily observed (Athens et al., [4]). It then seems rather unlikely that the effect would be of a transient nature.
With this in mind, we randomly assign the rather large half-life of a week (168 h), and define a prior with a CV of 200%. This corresponds to a very flat prior, and the information present in the experimental data will completely dominate the posterior. In our experience, the assignment of a completely uninformative prior, for example, a uniform prior ranging from zero to infinity would produce a very inefficient MCMC sampling and prolong run times considerably. In addition, when the disease state in the present model is taken into account, a half time of one week seems perfectly reasonable. If no recurrent strokes occur, most patients stabilize within a week and show no significant clinical changes beyond three months after stroke occurrence. With this in mind, a half-life for the downregulation of the increase in neutrophil plasma residence within the suggested approximate 95% confidence interval of 53–510 h is likely.
5.17. Power Term of Stem Cell Production (δ)
A semimechanistic pharmacodynamic model for neutrophil progression in conjunction with chemotherapy-induced myelosuppression has been developed (Friberg et al., [24]). In that model, a similar power term was employed for the feedback signal in the regulation of stem cell productions. The power term was different for each of six drugs being studied. The mean of these power terms is 0.172, with an SD of 0.0389. This corresponds to a CV of 0.22. Since this value is not based on any direct measurements of the neutrophil feedback signal, we double the uncertainty to 0.44.
REFERENCES
1 Alexanian, R., Donohue, D. M. (1965). Neutrophilic granulocyte kinetics in normal man. J. Appl. Physiol.70:803–808.[CSA]
2 Athens, J. W., Haab, O. P., Raab, S. O., Mauer, A. M., Ashenbrucker, H., Cartwright, G. E., Wintrobe, M. M. (1961a). Leukokinetic studies. IV. The total blood, circulating and marginal granulocyte pools and the granulocyte turnover rate in normal subjects. J. Clin. Invest.40:989–995.[INFOTRIEVE][CSA]
3 Athens, J. W., Raab, S. O., Haab, O. P., Mauer, A. M., Ashenbrucker, H., Cartwright, G. E., Wintrobe, M. M. (1961b). Leukokinetic studies. III. The distribution of granulocytes in the blood of normal subjects. J. Clin. Invest.40:159–164.[INFOTRIEVE][CSA]
4 Athens, J. W., Haab, O. P., Raab, S. O., Boggs, D. R., Ashenbrucker, H., Cartwright, G. E., Wintrobe, M. M. (1965). Leukokinetic studies. XI. Blood granulocyte kinetics in polycythemia vera, infection, and myelofibrosis. J. Clin. Invest.44:778–788.[INFOTRIEVE][CSA]
5 Bishop, C. R., Athens, J. W., Boggs, D. R., Warner, H. R., Cartwright, G. E., Wintrobe, M. M. (1968). Leukokinetic studies. XIII. A non-steady-state kinetic evaluation of the mechanisms of cortisone-induced granulocytosis. J. Clin. Invest.47:249–260.[INFOTRIEVE][CSA]
6 Bois, F. Y. (2000a). Statistical analysis of Clewell et al. PBPK model of trichloroethylene kinetics. Environ. Health Persp. Suppl.108:307–316.[CSA]
7 Bois, F. Y. (2000b). Statistical analysis of Fisher et al. PBPK model of trichloroethylene kinetics. Environ. Health Persp. Suppl.108:275–282.[CSA]
8 Bois, F. Y. (2001). Applications of population approaches in toxicology. Toxicol Lett.120:385–394.[INFOTRIEVE][CSA][CROSSREF]
9 Bois, F. Y., Maszle, D. R. (1997). MCSim: A Monte Carlo simulation program. J. Stat. Software [electronic publication]2:1–60.[CSA]
Bois, F. Y., Gelman, A., Jiang, J., Maszle, D. R., Zeise, L., Alexeef, G. (1996a). Population toxicokinetics of tetrachloroethylene. Arch. Toxicol.70:347–355.[INFOTRIEVE][CSA][CROSSREF]
Bois, F. Y., Jackson, E. T., Pekari, K., Smith, M. T. (1996b). Population toxicokinetics of benzene. Environ. Health Persp.104:1405–1411.[CSA]
Cartwright, G. E., Athens, J. W., Wintrobe, M. M. (1964). The kinetics of granulocytopoiesis in normal man. Blood24:780–803.[INFOTRIEVE][CSA]
Chaloner, K. (1996). Elicitation of prior distributions. In: Berry, D. A., Stangl, D. K., eds. Bayesian Biostatistics. New York: Marcel Dekker, pp. 141–156.
Cronkite, E. P., Fliedner, T. M. (1964). Granulocytopoiesis. N Engl. J. Med.270:1347–1352.[INFOTRIEVE][CSA]
Cronkite, E. P., Vincent, P. C. (1969). Granulocytopoiesis. Ser. Haemat.2:3–43.[CSA]
Cronkite, E. P., Burlington, H., Chanana, A. D., Joel, D. D., Reincke, U., Stevens, J. (1977). Concepts and observations on the regulation of granulocyte production. In: Baum, S. J., Ledney, G. D., eds. Experimental Hematology Today. New York: Springer-Verlag, pp. 41–49.
Dancey, J. T., Deubelbeiss, K. A., Harker, L. A., Finch, C. A. (1976). Neutrophil kinetics in man. J. Clin. Invest.58:705–715.[INFOTRIEVE][CSA]
Deubelbeiss, K. A., Dancey, J. T., Harker, L. A., Finch, C. A. (1975). Marrow erythroid and neutrophil cellularity in the dog. J. Clin. Invest.55:825.[CSA]
Donohue, D. M., Reiff, R. H., Hansen, M. L., Belson, Y., Finch, C. A. (1958). Quantitative measurements of the erythrocytic and granulocytic cells of the marrow and blood. J. Clin. Invest.37:1571–1576.[INFOTRIEVE][CSA]
Dresch, C., Najean, Y., Bauchet, J. (1975). Kinetic studies of 51Cr and DF32P labeled granulocytes. Br. J. Haematol.29:67–80.[INFOTRIEVE][CSA]
Ferber, G. (2002). Biomarkers and proof of concept. Methods Find Exp. Clin. Pharmacol.24Suppl C:35–40.[INFOTRIEVE][CSA]
Finch, C. A., Harker, L. A., Cook, J. D. (1977). Kinetics of the formed elements of human blood. Blood50:699–707.[INFOTRIEVE][CSA]
Fliedner, T. M., Cronkite, E. P., Killman, S. A., Bond, V. P. (1964). Granulocytopoiesis. II. Emergence and pattern of labeling of neutrophilic granulocytes in humans. Blood24:683–700.[INFOTRIEVE][CSA]
Friberg, L. E., Henningsson, A., Maas, H., Nguyen, L., Karlsson, M. O. (2002). Model of chemotherapy-induced myelosuppression with parameter consistency across drugs. Journal of Clinical Oncology20:4713–4721.[INFOTRIEVE][CSA][CROSSREF]
Galbraith, P. R., Valberg, L. S., Brown, M. (1965). Patterns of granulocyte kinetics in health, infection and in carcinoma. Blood25:683–692.[INFOTRIEVE][CSA]
Gelman, A. (2006). Prior distributions for variance parameters in hierarchical models. Bayesian AnalysisAccepted for publication. Manuscript available as pdf at http://www.stat.columbia.edu/∼gelman/research/published/taumain.pdf[CSA]
Gelman, A., Rubin, D. (1992). Inference from iterative simulation using multiple sequences (with discussion). Stat. Sci.7:457–511.[CSA]
Jonsson, E. N., McIntyre, F., James, I., Krams, M., Marshall, S. (2005). Bridging the pharmacokinetics and pharmacodynamics of UK-276,279 across healthy volunteers and stroke patients using a mechanistically based model for target mediated disposition. Pharm. Res.22:1236–1246.[INFOTRIEVE][CSA][CROSSREF]
Jonsson, F., Johanson, G. (2001a). A Bayesian analysis of the influence of GSTT1 polymorphism on the cancer risk estimate for dichloromethane. Toxicol. Appl. Pharmacol.174:99–112.[INFOTRIEVE][CSA][CROSSREF]
Jonsson, F., Johanson, G. (2001b). Bayesian estimation of variability in adipose tissue blood flow in man by physiologically based pharmacokinetic modeling of inhalation exposure to toluene. Toxicology157:177–193.[INFOTRIEVE][CSA][CROSSREF]
Jonsson, F., Johanson, G. (2002). Physiologically based modeling of the inhalation kinetics of styrene in humans using a Bayesian population approach. Toxicol. Appl. Pharmacol.179:35–49.[INFOTRIEVE][CSA][CROSSREF]
Jonsson, F., Bois, F. Y., Johanson, G. (2001). Physiologically based pharmacokinetic modeling of inhalation exposure of humans to dichloromethane during moderate to heavy exercise. Toxicol. Sci.59:209–218.[INFOTRIEVE][CSA][CROSSREF]
Lala, P. K., Maloney, M. A., Patt, H. M. (1964). A comparison of two markers of cell proliferation in bone marrow. Acta Haemat.31:1–8.[INFOTRIEVE][CSA]
Lunn, D. J., Aarons, L. (1998). The pharmacokinetics of saquinavir: a Markov chain Monte Carlo analysis. J. Pharmacokin. Biopharm.26:47–74.[CSA]
Lunn, D. J., Best, N., Thomas, A., Wakefield, J., Spiegelhalter, D. J. (2002). Bayesian analysis of population PK/PD models: general concepts and software. J Pharmacokinet Pharmacodyn.29:271–307.[INFOTRIEVE][CSA][CROSSREF]
Mary, J. Y. (1984). Normal human granulopoiesis revisited. II. Bone marrow data. Biomedicine & Pharmacotherapy38:66–77.[INFOTRIEVE][CSA]
Mauer, A. M., Athens, J. W., Ashenbrucker, H., Cartwright, G. E., Wintrobe, M. M. (1960). Leukokinetic studies. III. A method for labeling granulocytes in vitro with radioactive diisopropylfluorophosphate (DFP32). J. Clin. Invest.39:1481–1486.[INFOTRIEVE][CSA]
McMillan, R., Scott, J. L. (1968). Leukocyte labeling with 51chromium I. Technic and results in normal subjects. Blood32:738–754.[INFOTRIEVE][CSA]
Olsson Gisleskog, P., Karlsson, M. O., Beal, S. L. (2002). Use of prior information to stabilize a population data analysis. J. Pharmacokinet. Pharmacodyn.29:473–505.[CSA][CROSSREF]
Osgood, E. E., Seaman, A. J. (1944). The cellular composition of normal bone marrow as obtained by sternal puncture. Physiol. Rev.24:46–69.[CSA]
Sheiner, L. B., Wakefield, J. (1999). Population modelling in drug development. Stat. Methods Med. Res.8:183–193.[INFOTRIEVE][CSA][CROSSREF]
Smith, A. F. M. (1991). Bayesian computational methods. Phil. Trans. R. Soc. Lond.A:369–386.[CSA]
Smith, B. J. (2003). Bayesian Output Analysis Program (BOA), Version 1.0.1 for S-Plus and R.
Spear, R. C., Bois, F. Y. (1994). Parameter variability and the interpretation of physiologically based pharmacokinetic modeling results. Environ. Health Persp.102:61–66.[CSA]
Spiegelhalter, D. J., Abrams, K. R., Myles, J. P. (2004). Bayesian Approaches to Clinical Trials and Health-Care Evaluation. Chichester: John Wiley & Sons.
Uchida, T. (1971). Leukokinetic studies in peripheral blood. I. Neutrophilic granulocyte kinetics in normal man. Acta Haemat. Jap.34:164–185.[INFOTRIEVE][CSA]
Walker, R. I., Willemze, R. (1980). Neutrophil kinetics and the regulation of granulopoiesis. Reviews of Infectious Diseases2:282–292.[INFOTRIEVE][CSA]
Vaughan, S. L., Brockmyre, F. (1947). Normal bone marrow as obtained by sternal puncture. Blood1S:54–59.[CSA]
Footnotes
<sups>a</sups>For this parameter, no prior information existed and it was either fixed to zero ("informative" prior), given a value corresponding to a lognormal population variability of 30%, with an uncertainty of 1% ("semi-informative" prior), or assigned a value corresponding to a variability of 30% with an uncertainty of 200% ("uniformative" prior).
<sups>a</sups>All estimates based on the pooled output from the last 25 000 iterations of five Markov chains at convergence.<sups>b</sups>For this parameter, no prior information existed and it was either fixed to zero ("informative" prior), given a value corresponding to a lognormal population variability of 30%, with an uncertainty of 1% ("semi-informative" prior), or assigned a value corresponding to a variability of 30% with an uncertainty of 200% ("uninformative" prior).
<sups>a</sups>For this parameter, no prior information existed and it was fixed to zero.
By Fredrik Jonsson; E.Niclas Jonsson; FrédéricY. Bois and Scott Marshall
Reported by Author; Author; Author; Author