Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Pharmacokinet Pharmacodyn. Author manuscript; available in PMC 2010 June 8.
Published in final edited form as:
PMCID: PMC2882108

Pharmacodynamic analysis of time-variant cellular disposition: reticulocyte disposition changes in phlebotomized sheep


Most pharmacodynamic (PD) models of cellular response assume a time-invariant (i.e., constant) cellular disposition despite known changes in the disposition with time, such as the reticulocyte residence time in the systemic circulation during stress erythropoiesis. To account for changes in cellular disposition, a comprehensive PD model that involves endogenous erythropoietin (Epo), reticulocytes, and hemoglobin responses was developed in phlebotomized sheep that considers a time-variant reticulocyte residence time and allows for the simultaneous determination of changes in the cellular disposition and cellular production. Five sheep were phlebotomized to hemoglobin concentrations of approximately 4 g/dl. Epo concentrations, reticulocytes, and hemoglobin concentrations were frequently sampled for 5–7 days prior to and 25–30 days following the phlebotomy. Initial steady-state conditions were assumed and the time-variant reticulocyte residence time in the systemic circulation was semiparametrically represented using a constrained spline function. Hemoglobin production was modeled using a Hill function via an effect site compartment. The initial steady state reticulocyte residence time in the systemic circulation was estimated as 0.477 (0.100) (mean (SD)) days, which maximally increased 2.01- to 2.64-fold higher than the initial steady-state residence time 5.95 (0.899) days post-phlebotomy (P < 0.01). On average, the residence time returned to steady-state values 15.4 (2.36) days post-phlebotomy, which was not significantly different from the initial steady-state value (P > 0.05). The baseline hemoglobin production rate was estimated at 0.0929 (0.0472) g/kg/day and the maximum production rate under stress phlebotomy was estimated at 0.504 (0.0422) g/kg/day. These data indicate that endogenously released Epo under acute anemic conditions can increase hemoglobin production approximately 5-fold. The determined increase in reticulocyte residence time produced under stress erythropoiesis is similar to the commonly reported 2- to 3-fold increase observed in human patients.

Keywords: Reticulocyte maturation, Erythropoiesis, Hemapoietic cells, Hemoglobin production, Anemia, Time varient kinetics, Erythropoietin, Cellular transformations, Endogenous pharmacodynamics


Reticulocytes are immature red blood cells (RBCs) that initially reside in the bone marrow and subsequently are released into the systemic circulation where they develop into mature RBCs [1]. Reticulocytes are produced from erythroid progenitor cells (BFU-E and CFU-E) located primarily in the bone marrow. They are commonly used to monitor the “real-time” functional state of the erythroid bone marrow and play an important role in diagnoses and management of anemia [27]. Measurement of the reticulocyte count is valuable for determining if a patient has a functionally normal response to either endogenously produced or exogenously administered erythropoietin (Epo) [5,8]. Reticulocyte counts are also used to monitor bone marrow suppression by chemotherapy [9] and bone marrow engraftment following bone marrow transplantation [1013]. The maturation of erythroid progenitor cells into reticulocytes and ultimately RBCs is primarily controlled by Epo, a-35kD glycoprotein hormone produced by the pertibular cells of the kidney in response to oxygen need [1]. During the development from erythroid progenitor cells their hemoglobin content increases until it develops into a reticulocyte upon nucleus extrusion, where further maturation primarily involves the removal of ribosomal RNA, remodeling of the plasma membrane, and a progressive decrease in cell size, increasing the concentration of hemoglobin on a mass/volume basis, but not increasing the overall amount of hemoglobin within the cell (approximately 30 picograms per cell in adult humans) [1,3,4]. Therefore, a reticulocyte cell has approximately the same amount of hemoglobin as a RBC [14].

Generally, under normal erythropoietic conditions in humans (i.e., non-anemic or non-erythropoietically stimulated), reticulocytes reside in the systemic circulation approximately 24 h before maturation into RBCs [3,4]. However, during stress erythropoiesis (i.e., stimulated erythropoietic conditions), the reticulocyte residence time in the circulation increases to an estimated 2–3 days [2]. Reticulocytes produced under stress erythropoiesis also contain more residual ribosomal RNA, are larger, and have less flexible plasma membranes than those produced under normal conditions, and therefore are thought to be immature reticulocytes that under “normal” physiological conditions reside in the bone marrow until being released as more mature reticulocytes [2,4,14,15]. Accordingly, the reticulocyte counts increase under stress erythropoiesis not only due to increased reticulocyte production in response to Epo stimulation, but also due to a longer residence time in the systemic circulation. To account for the effect of this increased reticulocyte residence time in the systemic circulation on the reticulocyte count in clinical practice, one “rule of thumb” approach has been to divide the observed reticulocyte count or percentage by two to obtain a “more accurate” index of reticulocyte production when they are produced under stress erythropoiesis [4]. More sophisticated quantitative approaches utilize the maturity of the produced reticulocytes, referred to as the reticulocyte maturity index (RMI). The RMI is based on of individual reticulocyte RNA content commonly determined by fluorescence intensity measurements. However, the problem with existing methods utilizing RMI is that there is no standard method that corrects for the observed reticulocyte counts based on the RMI. Furthermore, these methods suffer from the difficulty in standardizing and cross-correlating across different instruments and laboratories [3].

Although different terminology is often used to describe reticulocyte or cellular residence times, in this manuscript it is defined as the time a cell resides in the sampling pool (the systemic blood circulation). Analogous terminology used by other authors includes “lifespan” and “maturation time” [1623]. Previous Epo-reticulocyte pharmacokinetic/pharmacodynamic (PK/PD) models have commonly assumed that all cells have a time invariant and common lifespan or maturation time (i.e., “point distribution”) [1621]. Recently, Epo-reticulocyte PK/PD models have been presented that allow for a fixed (time invariant) distribution of cell lifespans [22]. The time invariant lifespan assumption has been partially addressed by allowing for two different reticulocyte maturation times, a baseline and a stress-erythropoiesis maturation [23]. Other models that partially address the stress erythropoiesis effects on increasing reticulocyte lifespan incorporate an erythropoietin dose dependent lifespan [22]. However, the dynamic changes in reticulocyte residence time in the systemic circulation throughout the time course of stress erythropoiesis were not determined in any of the above studies. Given the known change in reticulocyte residence time in the systemic circulation following stress erythropoiesis, the clinical importance of the reticulocyte count/residence time change, and the limited attempts to quantify these changes; the primary objective of the present analysis was to dynamically determine the changes in reticulocyte residence time due to phlebotomy induced stress erythropoiesis. To achieve this aim, a time-variant cellular disposition model was developed and integrated into a comprehensive pharmacodynamic analysis of Epo, reticulocyte, and hemoglobin response in phlebotomized sheep.

Theoretical considerations

Unit response function

The unit response (UR) function of a linear cellular disposition, (Fig. 1, panel A), for hematopoietic cells may be described as:


where ‘a’ is the time it takes for an individual cell to appear in the systemic circulation (the sampling compartment) following progenitor stimulation, ‘b’ is the time the subsequent cell exits, and ‘t’ is measured, as ‘a’ and ‘b’, relative to the time of stimulation (time ‘s’ in Fig. 1, panel A). ‘U’ is the unit step function defined by:


Fig. 1
Unit Response, UR(t) (Eq. 1, Eq. 14), characterizing the disposition function of a hematopoietic cell. The cell is present in the systemic circulation when UR(t) = 1. Thus, ‘a’ and ‘b’ represent the disposition function ...

Therefore, UR(t) takes on the value of 1 when the cell is present, and a value of 0 when the cell is not present in the sampling compartment.

Time invariant cellular disposition

Consider hematopoietic progenitor cells, which in the following will be referred to simply as “cells”, that are to be stimulated over time by a stimulation function, fstim(t), which typically depends on endogenous growth factor(s) and/or exogenous drug. The dependence of fstim on these factors will be considered subsequently. The units for fstim are numbers of cells stimulated per unit time. Let it be assumed that these cells act fairly independently of each other, then one can apply the superposition principle. In the case where all cells have the same time invariant disposition parameters (a and b, Fig. 1, panel A) it follows that the number of cells present in the systemic circulation, N(t), is given by:


where * denotes the convolution operator. The integration limit of −∞ in Eq. (3) is to be interpreted to consider “all prior history” of the system. The above equation assumes a point distribution of cell residence times (i.e., zero variance). In order to readily fit Eq. (3) to observed data, it is convenient to transform it into a differential equation form. Equation (3) can be separated into two components as:


where t0 is the initial observation time. If we define:



Then Eq. (4) becomes:


If steady-state (SS) conditions are assumed to exist prior to t0 so that:


where fstimSS is the steady-state stimulation rate, then:


Differentiating Eq. (9) with respect to time gives:


Evaluation and differentiation of Eq. (6), Appendix A, results in:


Differentiation of Eq. (7) and substitution of Eq. (10) and Eq. (11) results in:




Extension to time variant cellular disposition parameters

Next consider the case where the cellular disposition parameters are time variant. That is, the residence time in the systemic circulation, ‘ba’, is a function of time. For simplicity we consider the time variance of the residence time to be due only to the ‘a’ parameter (Fig. 1, panel B). In this case the UR function becomes:


Here ‘s’ is the time that the cell is stimulated (Fig. 1, panel B). The value of ‘a’ is considered variable and determined at the time the cell is stimulated and remains fixed at that value thereafter. Again, a point distribution of cell residence times, at a given stimulation time, is assumed.

For this time-variant case, Eq. (4) becomes:


If similar to before we assume steady-state conditions for tt0, then:


where aSS is the steady-state ‘a’ value. Due to the steady-state assumption the first term, Ig(t), becomes nearly identical to Eq. (5), and similarly to before we get:


However, this time we write:


For convenience in derivation, we define:




Due to the unit step function, U(t), the integrand of Eq. (20) will be zero in the integration region u = t0 to u = t for ‘u’ values for which tua (u) < 0. Thus, the switch between zero and non-zero integrand values will occur at u = x value(s) where ‘x’, that is a function of time, is the solution(s) to the equation:


Potentially, Eq. (22) may have multiple solutions for ‘x’ at the current time ‘t’. However, as evident from the graphical representation of the solution(s) to Eq. (22) given in Fig. 2, if the a(t) never has a slope, a′(t), less than or equal to −1 then Eq. (22) only has one solution (Fig. 2, time ‘t1’). For a more detailed proof, see Appendix B. Accordingly, Eq. (20) becomes:



Fig. 2
Illustration of the relationship between ‘t’, ‘x(t)’, ‘a(t)’ and the condition for multiple solutions to Eq. 22. The scale on the ‘a(t)’ and the time axes are the same. The slope of the dashed ...

It is recognized from Eq. (22) that ‘x’ is simply the time for the progenitor cell stimulation of those “new” cells that are currently entering the systemic circulation at time ‘t’.

Differentiation of Eqs. (23) and (24) gives:



Eqs. (25) and (26) generalize to:


Similar to Eq. (6), evaluation and differentiation of Eq. (21) generalizes to:


Differentiation of Eq. (19) and substitution of Eqs. (27) and (28), followed by subsequent substitution of the resulting equation along with Eq. (17) into the time-variant equivalent of differentiated Eq. (7) results in:







Correction for cells removed by phlebotomy

Let FR denote the fraction of cells removed momentarily after the time of the phlebotomy, denoted ‘T’, then the equation to be integrated changes to a new “initial” condition at time ‘T’ according to:


where ε denotes an infinitesimally small time increment. However, just a change in a boundary condition according to Eq. (34) is insufficient to correct for the cells mechanically removed. The differential equation, Eq. (29), being integrated across the phlebotomy boundary point (t = T) must also be changed. If the original equation (Eq. (29)) is not changed then this equation will for some time period beginning immediately after the phlebotomy underestimate the number of cells because it contains an elimination component (i.e., Ih2(t)) that includes those cells mechanically removed at time ‘T’, which should not be included in the elimination. To correct for this, one needs to add a term that accounts for the elimination rate of those cells removed at time ‘T’ identical to their elimination had they not been removed by phlebotomy. Accordingly, the correction term (Appendix C) to be added is:


In deriving Eq. (35), it is recognized that the oldest cells in the removed cell population would be eliminated in vivo at t = T with a (projected) rate of FR · fstim (Tb), while the youngest cells would be eliminated at t = T + ba(x(T)) with a rate of FR · fstim (Ta(x(T))). Thus, the correction term is active from the time period t = T to t = T + ba(x(T)), consistent with the above considerations.

Summary of general key equation

For optimal readability the general key equation, which follows from Eqs. (29) and (35), is produced below together with support equations, although the latter equations are reproduced from above:



x(t)=ta(x(t))x(t)=11+a(x(t)),t>t0 and a(t)>1N(t)=N(t0)=fstimSS·(baSS),tt0x(t)=taSS,tt0a(t)=a(t0)=aSS,tt0fact(t)=fact(t0)=fstimSS,tt0b>a(t)0N(t)={N(t)fort=T(1FR)·N(T)fort=T+εU(t)={1ift00otherwise

Specific modeling of time variant disposition

The above Eq. (36) describes the general case considering any arbitrary variation over time in the a(t) “parameter”. The following section is aimed at specifically modeling the a(t) variation over time.

The time-variant cellular disposition parameters are not readily determined by modeling cellular count data from the cells of interest alone. The observed response is a function of both the stimulation rate and the disposition “parameters”, which are both changing as functions of time, causing difficulty in the estimation procedure. One possible solution is to utilize shared disposition parameters or functions from “coupled” cell data (i.e., cells that transform into each other) to assist in the simultaneous determination of the stimulation rate and the disposition “parameters”.

To determine the time-variant reticulocyte residence time, the present analysis utilizes the shared input and disposition functions of reticulocytes and hemoglobin to overcome the above problem, since a reticulocyte has approximately the same amount of hemoglobin as a mature RBC [1]. Accordingly, both the reticulocyte count and hemoglobin response share the same lag time, a(t), for appearance in the systemic circulation. The release of the hemoglobin into the systemic circulation is assumed to enter the circulation within a reticulocyte before it matures into an RBC. The ‘b’ parameter of the reticulocyte disposition was assumed to be time invariant. Furthermore, both the reticulocytes and hemoglobin share the stimulation function, fstim, which is multiplied by a conversion constant, KR, to change from units of hemoglobin to number of reticulocytes (i.e., has units of 103 cells/g of Hb). Additionally, since the same fractional removal in cells will occur for RBCs (i.e., hemoglobin) and reticulocytes, both share the FR parameter as well. Thus, reticulocytes and hemoglobin are “coupled” through several shared parameters and functions which generally assist in the analysis.

The stimulation function

Measurement of hematopoietic cells in vivo is done in terms of concentrations, such as number of erythroid cells or grams of hemoglobin per volume of blood, therefore the number of cells (or grams of hemoglobin), ‘N’, is divided by the volume, ‘V’. However, in the analysis, due to the measurements units, the ‘V’ fuses with the parameters of the stimulation function, fstim, so that it is not estimated directly in the regression analysis. Erythropoietin was considered the stimulator for the stimulation of the progenitor cells. Accordingly, Ce is the concentration of Epo at the biophase (effect site). The analysis indicated that Ce is kinetically distinct from the Epo plasma concentration, Cp, and the following convolution relationship:


was used to model the plasma/biophase relationship resulting in the following differential equation:




where CeSS and CpSS denote the steady-state effect site and plasma Epo concentrations, respectively, and CpSS was set equal to the initial (first observation) fitted plasma Epo concentration. The stimulation of the erythroid progenitor cells (quantified in terms of hemoglobin) was related to Ce by the Hill equation:




where qstim is the redefined stimulation function that depends on Ce, Emax is the maximal hemoglobin stimulation rate in g/dl/day, EC50 is the effect site Epo concentration that results in 50% of maximal hemoglobin stimulation rate (Emax), and qstim (CeSS) is the redefined steady-state stimulation rate.

Specific PD model

The analysis involved simultaneously fitting to both the hemoglobin and reticulocyte count concentration versus time profiles, denoted Hb(t) and R(t), respectively. Let ‘d’ denote the time the RBCs are removed from the systemic circulation relative to progenitor stimulation (i.e., the RBC equivalent ‘b’ parameter from Eq. (14) and Fig. 1) and ‘b’ denote the time the reticulocytes mature into RBCs relative to the stimulation. Both the ‘b’ and ‘d’ parameters are considered time invariant. As before, let a(t) denote the time that reticulocytes enter into the systemic circulation relative to progenitor stimulation. The time variance of a(t) is shared by the hemoglobin kinetics because the release of hemoglobin into the systemic circulation is assumed to initially enter within a reticulocyte. To ensure that b > a(t) ≥ 0 for all ‘t’ and for computational convenience, it is helpful to parameterize a(t) as:


Therefore, τ(t) is the (time-variant) reticulocyte residence time in the systemic circulation. The parameterization of Eq. (36) using Eq. (42) (Appendix D) with the appropriate substitutions (i.e., Eqs. (40) and (41)) from above results in the following specific key equations:










The final model is schematically detailed in Fig. 3.

Fig. 3
PD model schematic. Symbols are defined in the Glossary

Materials and methods


All animal care and experimental procedures were approved by the University of Iowa Institutional Animal Care and Use Committee. Five healthy young adult sheep ranging in age from 2- to 4-months old and weighing 21.2 (3.60) kg (mean (SD)) were selected. Animals were housed in an indoor, light- and temperature-controlled environment, with ad lib access to feed and water. Prior to study initiation, jugular venous catheters were aseptically placed under pentobarbital anesthesia. Intravenous ampicillin (1 g) was administered daily for 3 days following catheter placement.

Study protocol

Blood samples (~0.5 ml/sample) for plasma Epo, reticulocyte counts, and hemoglobin were collected for 5–7 days to determine baseline values prior to conducting a single controlled phlebotomy over several hours to induce acute anemia. Animals were phlebotomized to hemoglobin concentration of approximately 4g/dl. To maintain a constant blood volume during the procedure, equal volumes 0.9% NaCl solution were infused for each volume of blood removed. Blood samples were subsequently collected at least daily for 25–30 days post-phlebotomy. No iron supplementation other than that in the animal’s feed was given. Plasma iron and total iron binding capacity (TIBC) were monitored to ensure animals did not become iron deficient. To minimize hemoglobin and red cell loss due to frequent blood sampling, blood was centrifuged, the plasma removed, and the red cells re-infused.

Sample analysis

Plasma Epo concentrations were measured in triplicate using a double antibody radioimmunoassay (RIA) procedure as previously described (lower limit of quantitation 1mU/ml) [24]. All samples from the same animal were measured in the same assay to reduce variability. The reticulocyte count was determined by flow cytometry (FACScan, Bection–Dickinson, San Jose, CA) as previously described [25]. Hemoglobin concentrations were measured spectrophotometrically using an IL482 CO-oximeter (Instrumentation Laboratories, Watham, MA). In total, 30–45 samples were analyzed per variable per subject.

PD modeling details

The proposed PD model is schematically detailed in Fig. 3. The pharmacokinetics of Epo (i.e., Epo plasma concentrations) were non-parametrically represented using a generalized cross-validated cubic spline function, with nodes at every data point [26]. The hemoglobin and reticulocyte-plasma Epo concentration relationship was modeled by simultaneously fitting Eqs. (38)(40) and Eqs. (43)(50) to the modeled plasma Epo, observed blood hemoglobin, and observed reticulocyte count concentration–time data of each subject. The lag time parameter, ‘d’, between progenitor cell stimulation and RBC removal from the systemic circulation was fixed at ‘aSS + red cell (RBC + reticulocyte) lifespan’. The normal red cell lifespan has previously been determined in sheep using [14C] cyanate label and found to be 114 days [27]. A fifth order constrained spline was used to semiparametrically estimate the residence time of the reticulocytes, τ(t), as a function of time. The residence time by the spline function was constrained to remain at the pre-phlebotomy steady-state value (τ0 = τSS) until some time (T1) after the phlebotomy (τ′(t) = 0 for tT1 < T2) and then return to a steady-state value (τend) at some later point (T2, τ′(t) = 0 for tT2), as detailed in Appendix E.

Computational details

All modeling was conducted using WINFUNFIT, a Windows (Microsoft) version evolved from the general non-linear regression program FUNFIT [28], using weighted least squares. Data points were weighted by yobs2, where yobs is the observed hemoglobin or reticulocyte concentration. The analysis required the numerical solution to delay differential equations (DDE) because of the lag times involved in the mathematical model. This was done by WINFUNFIT using a DDE solver module that is based on a adaptation of the code “RETARD” [29] (see electronic supplementary material). The fraction of hemoglobin and reticulocytes removed, FR, was estimated by WINFUNFIT using an generalized events processing module that allows parameterized events to change during the iterative parameter estimation procedure. Accordingly, the magnitude and time for the events and delays can be estimated simultaneously with other unknown parameters of the mathematical model. Akaike’s Information Criterion (AIC) was used to compare the time variant model with 11 estimated parameters (four cellular production parameters (Emax, EC50, ke, and FR) and seven residence time parameters/coefficients (‘b’ and the spline parameters/coefficients)) to the identical time invariant model with only six estimated parameters (four cellular production and two residence time parameters (‘a’ and ‘b’)) [30].

To summarize the uncertainty in the individual subject parameter estimates, the mean percent standard error (MSE%) of the estimate was calculated for each parameter as:


where SEi and Pi are the standard error of the parameter and the estimate of the parameter for the ith subject, respectively, and ‘n’ is the number of subjects.

Statistical analysis

Mean reticulocyte residence times in the systemic circulation were compared using paired two tailed Student’s t-test at selected points in time using Microsoft® Excel 2002 SP3 (Microsoft Corporation, Redmond, WA). Statistically significant differences were determined at the α = 0.05 type I error rate.


The simultaneous fit to the Epo, reticulocyte, and hemoglobin data for two representative subjects is shown in Fig. 4. Good agreement between the observed and predicted values is evident from the fits (overall r = 0.978 (0.0084) (mean (SD)), hemoglobin r = 0.981 (0.0203), reticulocyte r = 0.966 (0.0114)). The PD model parameters are summarized in Table 1. The mean determined maximal rate of hemoglobin production, Emax, was estimated to be 1.29 (0.471) g/dl/day. The mean lag time between progenitor stimulation and reticulocyte disappearance from the systemic circulation, ‘b’, was estimated at 1.24 (0.167) days and the mean pre-phlebotomy steady-state reticulocyte residence time in the systemic circulation, τ0, was estimated to be 0.447 (0.100) days. Accordingly, the mean pre-phlebotomy steady-state lag time between progenitor stimulation and reticulocyte (hemoglobin) appearance in the systemic circulation, ‘aSS’, was estimated at 0.797 (0.193) days.

Fig. 4
Representative individual subject fits of model equations to observed plasma erythropoietin ([diamond with plus]), hemoglobin (□), and reticulocyte (+) count concentrations (panels A and B are different subjects).
Table 1
Parameter estimates for the time-variant reticulocyte residence time model (n = 5)

The reticulocyte systemic circulation residence time determined for each subject rapidly increased shortly after the phlebotomy, reaching a maximal value of 0.761–1.24 days at 5.95 (0.899) days post-phlebotomy (Fig. 5). Figure 5 displays the residence time for cells currently being stimulated, not the residence of reticulocytes currently entering the bloodstream or circulating in the systemic blood. The time post-phlebotomy that the residence time began to deviate from steady-state conditions was estimated at 0.730 (0.834) days. The maximal residence time was 2.01–2.64-fold higher than the initial steady-state residence time (P < 0.01). Following establishment of the maximal residence time, the residence time rapidly dropped in all subjects below its pre-phlebotomy value, to a minimum value of 0.471–0.906-fold of the initial steady-state residence time (P < 0.05). It then remained at that value or increased to a new steady-state value. On average, the residence time returned to steady-state conditions 15.4 (2.36) days following the phlebotomy. In all but one subject, the ending steady-state reticulocyte residence time, τend, was near the initial value (P > 0.25, including and excluding the single subject with the ending residence time well above the initial value). The reason for this different behavior in the single subject is not clear, but it is evident from the subject’s data (Fig. 4, panel B) that the reticulocyte count was substantially elevated compared to the baseline count well past the phlebotomy, supporting that the reticulocyte residence time likely remained elevated.

Fig. 5
Determined reticulocyte residence times in the systemic circulation for the individual subject fits to the model equations (n = 5). The residence times were represented using a semiparametric spline function and are displayed as the residence times, τ( ...

The MSE% for the primary parameter estimates are also displayed in Table I. The other residence time parameters (see Appendix E), T1 and T2, had MSE% of 1.2% and 1.4%, respectively. The two unconstrained coefficients of the residence time fifth order spline function, α4 and α5, had MSE% of 0.7% and greater than 1000%, respectively. The high MSE% for α5 was due to two subjects, as the other three subjects had a MSE% of 0.6%. Modeling the data instead with a fourth-order residence time spline function still resulted in a good model fit to the data for the two subjects with the high SE for the α5 coefficient, however, it resulted in an unacceptably poor fit to the data from the other three subjects. Therefore, a fifth-order spline function was used to model the time-variant residence time as the minimal model that resulted in an acceptable fit to the all the subject data.

In all subjects, the time variant reticulocyte residence model was superior to the identical time invariant reticulocyte residence model, where a single residence time in the systemic circulation was estimated, based on AIC. The difference in the fits between these two models for a representative subject is displayed in Fig. 6.

Fig. 6
Representative subject fit to observed hemoglobin (□) and reticulocyte (+) count concentrations using a time variant or time invariant reticulocyte residence time model (panel A). The lines represent the predicted hemoglobin and the predicted ...


Cellular production

The determined steady-state lag time between effect-site Epo stimulation and reticulocyte or hemoglobin appearance in the systemic circulation, the ‘aSS’ value, of 0.797 (0.193) days is similar to the corresponding parameter we have determined in previous sheep phlebotomy studies of 0.72 days [23]. The current model includes an effect site compartment, while the previous analysis did not. In contrast to the previous analysis [23], this analysis also includes hemoglobin data in the estimation procedure. The determined ‘aSS’ value (0.797 days or 19.1 h) is approximately the same as a similar physiological parameter determined in cynomolgus monkeys of 14.95 h and in humans of 10.76 h [19, 20].

The determined CpSS value and qstim function for each subject estimates a steady-state hemoglobin production of 0.0929 (0.0472) g/kg/day, assuming a blood volume of 0.75 dl per kg bodyweight. Furthermore, the maximum production rate of hemoglobin following the observed Epo stimulation (i.e., maximum Ce) was estimated at 0.504 (0.0422) g/kg/day, indicating that endogenous Epo released under maximally stimulated acute anemic conditions can increase hemoglobin (RBC) production approximately 5-fold. The determined 5-fold increase in hemoglobin production following acute anemia in sheep is similar to the 3- to 5-fold increase reported in human patients [31].

Reticulocyte residence time

As would be predicted from the current understanding in the scientific community of reticulocyte maturation under stress erythropoiesis, the reticulocyte residence time in the systemic circulation increased post-phlebotomy (Figs. 5 and and6,6, panel B). The mean pre-phlebotomy steady-state reticulocyte residence time, τ0, was estimated at 0.447 (0.100) days (10.7 h) (Table 1). The estimated τ0 is similar to other estimates in sheep of 0.61 days and the commonly reported value of 1 day in humans [3, 23, 32]. From the steady-state value, the reticulocyte residence time increased 2.01- to 2.64-fold during stress erythropoiesis, before dropping below pre-phlebotomy values. The determined maximally increased reticulocyte residence time following a phlebotomy is approximately the same as that previously estimated in vivo in sheep and humans of 4.5-fold and 2–3-fold, respectively [23, 32]. It is also very close to the commonly recognized increase in reticulocyte residence time in the systemic in humans due to stress erythropoiesis of 2–3 days [2]. The drop of the mean reticulocyte residence time below the pre-phlebotomy steady-state value is not generally noted in the literature, and may be due to an overshoot/rebound phenomenon often observed in biological systems.

It is important to note that the determined reticulocyte residence time is not the reticulocyte lifespan. The latter is not possible to determine accurately in vivo by only sampling in the systemic circulation. A reticulocyte initially normally exists in the bone marrow of sheep and humans [1, 3, 32], and therefore cannot be observed until it is released from the marrow into the systemic circulation. Accordingly, the determined drop below baseline reticulocyte residence time post-phlebotomy most likely represents the release of more mature reticulocytes from the bone marrow into the circulation through some rebound/overshoot phenomenon.

Model details

Utilizing a time-variant reticulocyte residence time model was superior to a time invariant reticulocyte residence time model for all subjects based on the AIC (Fig. 6). As can also be seen in Fig. 6, an adequate simultaneous fit to both the reticulocyte and hemoglobin data is not possible with a time invariant reticulocyte residence time model. The current model is also an improvement over previous Epo, reticulocyte, and hemoglobin integrated PK/PD models published by our laboratory in accounting for hemoglobin production [16], since previous models did not include the hemoglobin in reticulocytes as part of the total measured hemoglobin. The hemoglobin contained in the reticulocyte, as previously modeled, only added to the measured hemoglobin upon maturation of the reticulocyte into a RBC. Furthermore, previous models used the percent reticulocytes instead of the absolute reticulocyte count, the former which is not as easily interpreted physiologically. The current analysis also provides an exact correction for the effect of removing erythroid cells during the phlebotomy process.

Because the current model shared functions through the improved hemoglobin accounting, a time-variant reticulocyte residence time could be modeled, which is widely known to occur physiologically, but apparently has not been previously modeled. It may be physiologically relevant to include time variance in the ‘b’ parameter in addition to the ‘a’ parameter. However, the time-variant residence could not be modeled as such because a variable ‘b’ value in addition to a variable ‘a’ value would not be possible to determine in a reliable way. Since reticulocytes are released from the bone marrow at a younger age following stress erythropoiesis, modeling the time variance due to changes in the ‘a’ value accounts for the primary physiological features of reticulocytes produced under stress erythropoiesis. The ability to estimate a time-variant ‘a’ value is primarily due to the simultaneous fitting to the “coupled” hemoglobin and reticulocyte responses. Two key features of the “coupled” responses are the longer residence time in the systemic circulation of hemoglobin compared to reticulocytes and the larger pool of hemoglobin relative to reticulocytes. Accordingly, a substantial change in the ‘a’ value, and thus the hemoglobin and reticulocyte input, has a relatively large and rapid effect on the observed reticulocyte response, but only a relatively minor effect on the observed hemoglobin response. Alternatively, a model with a time invarianta’ parameter but with a time variantb’ parameter could be formulated and fitted to the data to estimate the reticulocyte residence time in the systemic circulation. However, a model formulated as such would not acknowledge the primary features of erythropoiesis physiology. The model also determined a time-variant hemoglobin residence time in the systemic circulation. However, we only report and discuss the results of the model in terms of the change in the reticulocyte residence time, τ, and the observed reticulocyte response due to the small magnitude of the percent change in the hemoglobin residence time (<1%), the fixing of the ‘d’ parameter, and the limited observation time relative to the residence time of hemoglobin.

Another important model assumption is the fixed ‘d’ value of ‘aSS + red cell lifespan’ for all cells, since it has been experimentally demonstrated in mice and rats that red cells produced under stress erythropoietic conditions have a shortened lifespan or residence time in the systemic circulation [33, 34]. However, the lifespan of the red cells would need to be reduced by greater than 70% (i.e., to less than 30 days) to have an effect, since this was the longest period of time of observation post-phlebotomy. A final simplifying assumption is that a single, calculated hemoglobin to reticulocyte count conversion constant, KR, was used for all subjects. The KR parameter acts only as a scalar though, so any changes in the magnitude of the reticulocyte residence time would still be observed even if the value for KR varied from subject-to-subject. Additionally, a constant KR would be inappropriate if the amount of hemoglobin per erythroid cell varied with time, as may occur if subjects became iron deficient. However, monitoring of iron status did not indicate that any subjects became iron deficient throughout the time course of the experiments.

The model, as currently described, does not account for a distribution(s) of cell residences, since it was assumed that all cells stimulated at a given time have the same “pre-determined” residence time. However, previous work has demonstrated little advantage in moving to a more complicated residence or lifespan distribution model compared to a point distribution model due to problems estimating the standard deviation of the cell lifespans [22]. The simultaneous determination of a cellular distribution(s) and changes in the individual ‘a’ and/or ‘b’ values would be mutually confounding and would create tremendous estimation problems.

The present analysis did not indicate that a′(t) at any time fell at or below a value of −1 (nor did τ′(t) fall at or above a corresponding value of +1). Thus, the problem associated with accounting for multiple solutions, possible according to Eq. (22), fortunately did not arise in this analysis (see Fig. 2, showing multiple solutions when a′(t) becomes < −1). Physiologically, what occurs when multiple solutions to Eq. (22) arise is that cells that were stimulated at a later point in time are entering the sampling compartment prior to cells that were stimulated earlier (i.e., “younger” cells are entering prior to “older” cells), resulting in multiple groups (solutions) of cells that were stimulated at different times entering the systemic circulation at exactly the same time. Additionally, mathematically Eq. (30) becomes undefined when a′(t) = −1, because there is an infinite number of solutions along the finite range where a′(t) continuously equals −1. If the analysis had indicated that a′(t) at any time fell at or below −1, we would have had to resort to much more complex numerical solutions to the already complex delay-type differential equations describing the proposed model.

The time-variant cellular disposition model presented in the theoretical section is applicable to modeling any cellular response where changes in cellular residences is expected and may be extended to incorporate unique features of the specific physiological system being modeled. Furthermore, the specific Epo, reticulocyte, and hemoglobin PK/PD model that was presented could readily be applied to humans, both adult and neonate, and other species. Future modeling work and more “complete” data may also allow for the relaxation of certain model assumptions (i.e., a point distribution of cell residences), developing an even more general PD model of cellular responses.


In summary, by sharing parameters and functions between the measured hemoglobin and reticulocyte count response, the reticulocyte residence time in the systemic circulation was accurately determined as it changed with time. The utilized methodology allows for the dynamic determination of time-variant cellular disposition or residence time without directly measuring the cellular residence times. The determined increase in reticulocyte residence time from reticulocytes produced under stress erythropoiesis and the determined increase in hemoglobin production is consistent with other experimental data and with the scientific community’s current understanding of erythropoiesis physiology. Furthermore, the determined reduction below baseline residence time post-phlebotomy may serve to generate future experiments for determining why or how physiologically the reduction in reticulocyte residence times occurs in vivo in sheep, and ultimately in humans.

Supplementary Material

Supplemental material


The authors thank Dr. Robert T. Cook and the personnel of the Iowa City VAMC Pathology and Laboratory Medicine Services for assistance in performing the flow cytometric measurements of the reticulocytes. This work is supported by the United States Public Health Service, National Institute of Health grants PO1 HL46925, and R21 GM57367, and by the Veterans Administration Medical Center, Iowa City, Iowa.


Lag-time between stimulation and appearance of a cell in the sampling-compartment
Akaike’s Information Criterion
Lag-time between stimulation and disappearance of a cell from the sampling compartment
Burst Forming Unit-Erythroid
Effect site erythropoietin concentration
Plasma erythropoietin concentration
Colony Forming Unit-Erythroid
Lag-time between stimulation of an erythroid progenitor and disappearance of a RBC
Delay Differential Equation
Effect site erythropoietin concentration where 50% of Emax occurs
Maximal hemoglobin stimulation rate
Infinitesimally small time increment
Stimulation rate of cells as a function of time
Fraction of cells removed from phlebotomy
Hemoglobin concentration
First-order rate constant of the effect site concentration
Hemoglobin to reticulocyte count conversion constant
Mean percent standard error
Number of cells
Stimulation rate of cells as a function of Ce
Reticulocyte count (concentration)
Red Blood Cell
Reticulocyte Maturity Index
Time of stimulation
Standard deviation
Standard error
Steady-State (super- or sub-scripted with other terms)
Initial observation time
Time of phlebotomy
Time post-phlebotomy that steady-state τ values are deviated from
Time post-phlebotomy that steady-state τ values are returned to (T2 > T1)
τ = ba
Reticulocyte residence time in the systemic circulation
τ0 = baSS
Pre-phlebotomy (initial) steady-state reticulocyte residence time in the systemic circulation
End steady-state reticulocyte residence time in the systemic circulation
Arbitrary integration variable
Unit step function
Unit Response function
Time of stimulation of cells currently entering the sampling compartment
z = ln(b − τ(x))
Logarithmically transformed a

Appendix A

Derivation of Eq. (11)

By evaluation of Eq. (6):

Ih(t)=0, fortt0<a

Ih(t)=t0tafstim(u)du, foratt0<b

Ih(t)=tbtafstim(u)du, fortt0b>a

Differentiation of Eq. (A1) through (A3) gives:

Ih(t)=0, fortt0<a

Ih(t)=fstim(ta), foratt0<b

Ih(t)=fstim(ta)fstim(tb), fortt0b>a

Equation (A4) through (A6) generalize to:


which proves Eq. (11).

Appendix B

Proof of unique solution to Eq. 22 when a′(t) > −1.

Equation 22 may be written as:


with y [equivalent] x(t). Treating ‘t’ as fixed the Jacobian determinant, |JG|, of Eq. (B1) is:


Thus by the Implicit Function Theorem [35], Eq. 22 has a unique solution if |JG| ≠ 0, i.e., if 1 + a′(y) ≠ 0, which is satisfied when a′(y) > −1 or when a′(y) < −1 for all ‘y’. The condition a′(t) > −1 for all ‘t’ values ensures that a′(y) > −1 for all ‘y’ values, thus ensuring a unique solution to Eq. 22.

Appendix C

Derivation of Eq. (35)

It is recognized that the phlebotomy at time t = T can only remove cells already present in the circulation. If we separate Eq. (15) such as:


and define:



It is apparent that the Ij(t) term gives all the cells that have already entered into the circulation at time, t = T, and thus all the cells that can be removed by the phlebotomy.

By evaluation of Eq. (C2):




Since the Ij(t) term represents the cells that can be removed by the phlebotomy, it is multiplied by (1 − FR) of Eq. (34) at time t = T + ε, resulting in:

N(t)=(1FR)·Ij(t)+Ik(t)=Ij(t)+Ik(t)FR·Ij(t),  fort>T

The first term on the right hand side of Eq. (C7), Ij(t) + Ik(t), is recognized as N(t) of Eq. (15), which if differentiated will result in Eq. (29). Therefore, −FR · Ij(t) represents the correction term. Since tT > 0 in Eq. (C7), Ij(t) only has non-zero values from time t = T to t = T + ba (x (T)) as given by Eq. (C5). Differentiation of Eq. (C5) results in:


Thus, the correction term generalizes to:


which proves Eq. (35).

Appendix D

Alternative parameterization of Eq. (36)

A useful alternative parameterization to ensure that the conditions of Eq. (42) hold and for computational convenience follows: from Eq. 22 define:


The exponential transformation ensures that a (x (t)) > 0 for any x(t) value. Rearrangement and differentiation give:



From Eqs. (33) and (42), it follows that:



Substituting Eq. (D2) through Eq. (D5) into Eq. (36) results in:





Additionally, substitution of Eq. (D2) into the right-hand side of Eq. (D1) results in:


Differentiation of Eq. (D5) and (D9), with substitution and rearrangement results in the following additional conditions to Eq. (D6):




The above τ′(t) < 1 constraint corresponds to a′(t) > −1 to ensure that Eq. (22) has a single solution (Fig. 2 and Appendix B).

Appendix E

Residence time spline function

The reticulocyte systemic circulation residence time, τ(t), was modeled as follows:




Such that:







In order to satisfy the above constraints, the fifth-order polynomial that represents τ(t) when t1tt2 was parameterized in terms of the parametric terms, τ0, τend, T1 and T2, and the non-parametric coefficients, α4 and α5.

Contributor Information

Kevin J. Freise, Division of Pharmaceutics, College of Pharmacy, The University of Iowa, Iowa City, IA 52242, USA.

John A. Widness, Department of Pediatrics, College of Medicine, The University of Iowa, Iowa City, IA 52242, USA.

Robert L. Schmidt, Department of Pediatrics, College of Medicine, The University of Iowa, Iowa City, IA 52242, USA.

Peter Veng-Pedersen, Division of Pharmaceutics, College of Pharmacy, The University of Iowa, Iowa City, IA 52242, USA.


1. Hoffman R, Benz EJ, Jr, Shattil SJ, Furie B, Cohen HJ, Silberstein LE, McGlave P. Hematology: basic principles and applications. 4th edn. USA: Elsevier Inc; 2005.
2. Brugnara C. Use of reticulocyte cellular indices in the diagnosis and treatment of hematological disorders. Int J Clin Lab Res. 1998;28:1–11. [PubMed]
3. Brugnara C. Reticulocyte cellular indices: a new approach in the diagnosis of anemias and monitoring of erythropoietic function. Crit Rev Clin Lab Sci. 2000;37:93–130. [PubMed]
4. Hillman RS, Ault KA, Rinder HM. Hematology in clinical practice. 4th edn. USA: McGraw-Hill Companies, Inc; 2005.
5. Riley RS, Ben-Ezra JM, Tidwell A, Romagnoli G. Reticulocyte analysis by flow cytometry and other techniques. Hematol Oncol Clin North Am. 2002;16:373–420. vii. [PubMed]
6. Hirose A, Yamane T, Shibata H, Kamitani T, Hino M. Automated analyzer evaluation of reticulocytes in bone marrow and peripheral blood of hematologic disorders. Acta Haematol. 2005;114:141–145. [PubMed]
7. Watanabe K, Kawai Y, Takeuchi K, Shimizu N, Iri H, Ikeda Y, Houwen B. Reticulocyte maturity as an indicator for estimating qualitative abnormality of erythropoiesis. J Clin Pathol. 1994;47:736–739. [PMC free article] [PubMed]
8. Pradella M, Cavill I, d’Onofrio G. Assessing erythropoiesis and the effect of erythropoietin therapy in renal disease by reticulocyte counting. Clin Lab Haematol. 1996;18:35–37. [PubMed]
9. Aulesa C, Ortega JJ, Jou JM, Colomer E, Montero J, Olano E, Remacha AF, Martino R, Yoldi F, Escribano I. Flow cytometric reticulocyte quantification in the evaluation of hematologic recovery—Spanish multicentric study-group for hematopoietic recovery. Euro J Haematol. 1994;53:293–297.
10. Davis BH, Bigelow NC. Flow cytometric reticulocyte quantification using thiazole orange provides clinically useful reticulocyte maturity index. Arch Pathol Lab Med. 1989;113:684–689. [PubMed]
11. Batjer JD, Riddell K, Fritsma GA. Predicting bone-marrow transplant engraftment by automated flow cytometric reticulocyte analysis. Lab Med. 1994;25:22–26.
12. Dubner D, Perez MD, Barboza M, Sorrentino M, Robinson A, Gisone P. Prognosis and bone marrow recovery indicators in bone marrow transplantation after total body irradiation. Med Buenos Aires. 2002;62:555–561. [PubMed]
13. Noronha JFA, De Souza CA, Vigorito AC, Aranha FJP, Zulli R, Miranda ECM, Grotto HZW. Immature reticulocytes as an early predictor of engraftment in autologous and allogeneic bone marrow transplantation. Clin Lab Haematol. 2003;25:47–54. [PubMed]
14. Houwen B. Reticulocyte maturation. Blood Cells. 1992;18:167–186. [PubMed]
15. Jandl JH. Blood: textbook of hematology. 2nd edn. Little, Brown and Company, USA; 1996.
16. Veng-Pedersen P, Chapel, Schmidt RL, Al-Huniti NH, Cook RT, Widness JA. An integrated pharmacodynamic analysis of erythropoietin, reticulocyte, and hemoglobin responses in acute anemia. Pharm Res. 2002;19:1630–1635. [PubMed]
17. Chapel SH, Veng-Pedersen P, Schmidt RL, Widness JA. A pharmacodynamic analysis of erythropoietin-stimulated reticulocyte response in phlebotomized sheep. J Pharmacol Exp Ther. 2000;295:346–351. [PubMed]
18. Krzyzanski W, Jusko WJ, Wacholtz MC, Minton N, Cheung WK. Pharmacokinetic and pharmacodynamic modeling of recombinant human erythropoietin after multiple subcutaneous doses in healthy subjects. Eur J Pharm Sci. 2005;26:295–306. [PubMed]
19. Ramakrishnan R, Cheung WK, Farrell F, Joffee L, Jusko WJ. Pharmacokinetic and pharmacodynamic modeling of recombinant human erythropoietin after intravenous and subcutaneous dose administration in cynomolgus monkeys. J Pharmacol Exp Ther. 2003;306:324–331. [PubMed]
20. Ramakrishnan R, Cheung WK, Wacholtz MC, Minton N, Jusko WJ. Pharmacokinetic and pharmacodynamic modeling of recombinant human erythropoietin after single and multiple doses in healthy volunteers. J Clin Pharmacol. 2004;44:991–1002. [PubMed]
21. Krzyzanski W, Ramakrishnan R, Jusko WJ. Basic pharmacodynamic models for agents that alter production of natural cells. J Pharmacokinet Biopharm. 1999;27:467–489. [PubMed]
22. Krzyzanski W, Woo S, Jusko WJ. Pharmacodynamic models for agents that alter production of natural cells with various distributions of lifespans. J Pharmacokinet Pharmacodyn. 2006;33:125–166. [PubMed]
23. Al-Huniti NH, Widness JA, Schmidt RL, Veng-Pedersen P. Pharmacodynamic analysis of changes in reticulocyte subtype distribution in phlebotomy-induced stress erythropoiesis. J Pharmacokinet Pharmacodyn. 2005;32:359–376. [PubMed]
24. Widness JA, Veng-Pedersen P, Modi NB, Schmidt RL, Chestnut DH. Developmental differences in erythropoietin pharmacokinetics: Increased clearance and distribution in fetal and neonatal sheep. J Pharmacol Exp Therapeut. 1992;261:977–984. [PubMed]
25. Peters C, Georgieff MK, Alarcon PAd, Cook RT, Burmeister LF, Lowe LS, Widness JA. Effect of chronic erythropoietin administration on plasma iron in newborn lambs. Biol Neonate. 1996;70:218–228. [PubMed]
26. Hutchinson MF, deHoog FR. Smoothing noise data with spline functions. Numer Math. 1985;47:99–106.
27. Mock DM, Lankford GL, Burmeister LF, Strauss RG. Circulating red cell volume and red cell survival can be accurately determined in sheep using the [14C]cyanate label. Pediatr Res. 1997;41:916–921. [PubMed]
28. Veng-Pedersen P. Curve fitting and modelling in pharmacokinetics and some practical experiences with NONLIN and a new program FUNFIT. J Pharmacokinet Biopharm. 1977;5:513–531. [PubMed]
29. Hairer E, Norsett SP, Wanner G. Solving ordinary differential equations I: nonstiff problems. 2nd edn. Berlin: Springer-Verlag; 1993.
30. Akaike H. Automatic control: a new look at the statistical model identification. IEEE Trans. 1974;19:716–723.
31. Hillman RS, Finch CA. Erythropoiesis: normal and abnormal. Semin Hematol. 1967;4:327–336. [PubMed]
32. Hillman RS. Characteristics of marrow production and reticulocyte maturation in normal man in response to anemia. J Clin Invest. 1969;48:443–453. [PMC free article] [PubMed]
33. Shimada A. The maturation of reticulocytes. II. Life-span of red cells originating from stress reticulocytes. Acta Med Okayama. 1975;29:283–289. [PubMed]
34. Stohlman F., Jr Humoral regulation of erythropoiesis. VII. Shortened survival of erythrocytes produced by erythropoietine or severe anemia. Proc Soc Exp Biol Med. 1961;107:884–887. [PubMed]
35. Krantz SG, Parks HR. The implicit function theorem: history, theory, and applications. Boston: Birkhauser; 2002.