Subjects and hematological biomarker data

The analysis included data from a Phase 1, multicenter, randomized, double-blind, placebo-controlled, parallel-group clinical trial. The primary aim of this study was to investigate the hematological and visual effects, safety, and tolerability of bitopertin. A total of 67 healthy males and females aged < 50 years were enrolled in the trial. The subjects were administered bitopertin (10, 30, or 60 mg) or placebo orally once daily for 120 days, followed by a 120-day follow-up period.

Blood samples for hematological assessments were collected at baseline, week 1, week 2, every 2 weeks thereafter until week 16, week 17 (at treatment end), week 18, and every 2 weeks thereafter until the end of the follow-up period at week 34. Additionally, rich pharmacokinetic (PK) sampling and records of individual bitopertin doses were available. These data enabled the estimation of individual PK parameters using a population PK model (data not shown).

Development of the erythropoiesis modelPhysiological process of erythropoiesis

The initial stages of erythrocyte development occur in the bone marrow, where various precursor cells are found. These precursor cells undergo a series of maturation stages for approximately 5 days [10, 11]. The most mature precursor form is the reticulocyte. Reticulocytes reside in the bone marrow for approximately 3 days and in the blood for approximately 1 day, leading to a total lifespan of approximately 4 days [4, 10]. Hemoglobin synthesis occurs in the early reticulocyte maturation stages, mainly in the bone marrow. Hemoglobin synthesis is confined to these immature reticulocytes because, after this stage, they no longer possess the necessary cellular machinery and ribosomal RNA for protein synthesis [4]. The reticulocyte matures into an erythrocyte by losing its internal organelles and remodeling its plasma membrane. Under hematologically normal conditions, erythrocytes remain stay in circulation for approximately 120 days, before being engulfed by splenic macrophages [10, 12].

Integration of available biomarkers into the model development

Early-stage precursors in the bone marrow were not directly observed because of the invasive nature of the sampling. Thus, this work focused on observations from blood samples, particularly reticulocytes, and their subdivision into mature and immature forms. Both the reticulocyte total blood count (RET) and the fraction of immature reticulocytes (IRF) in the blood are indicative of erythropoietic activity in the bone marrow [4].

Hemoglobin is a critical biomarker of the blood’s oxygen transportation capacity. A hemoglobin deficiency is used as a clinical marker for anemia [13]. Hemoglobin is assessed through total hemoglobin blood concentration (Hbtot) and mean corpuscular hemoglobin (MCH) in erythrocytes.

Structural model

A previously developed model that focused on hemoglobin turnover was used as a framework upon which a more extensive erythropoiesis model could be built [9]. The model describes erythrocyte turnover and the mean hemoglobin content within the erythrocytes as two parallel chains of four transit compartments that all share the same parameter for the transit rate constant. Hemoglobin turnover was described as the net effect of several processes. First, there was the synthesis of hemoglobin inside the precursor cells. Second, there was the turnover of erythrocytes in the blood. Finally, a feedback mechanism increases precursor recruitment when Hbtot decreased. Bitopertin’s inhibitory effect on hemoglobin synthesis was modeled as an exposure-dependent decrease of the hemoglobin production rate. When erythrocytes with a reduced hemoglobin content enter the blood, Hbtot consequently decreases, inducing a homeostatic feedback that stimulated the precursor recruitment rate.

The previously described model was expanded based on the following assumptions: reticulocyte maturation pathways described using a four-compartment structure representing immature and mature reticulocytes in bone marrow and blood [14]; an equal reticulocyte maturation rate in bone marrow and blood [4]; and homeostasis (steady state) at baseline. Furthermore, to reflect swift reticulocyte dynamics, an empirical tolerance mechanism was evaluated [15].

The model was fitted simultaneously to the observed data on RET, erythrocyte count (RBC), MCH, and IRF. Individual AUCss values were used as the exposure metric driving the bitopertin drug effect [9].

Statistical model

The need to include inter-individual variability (IIV) on a model parameter was evaluated. The IIV in a parameter was modeled using a log-normal distribution \(^_}\). The random effect variable was assumed to follow a normal distribution \(_\sim N\left(0,_^\right)\) on parameter \(p\) for individual \(i\), where the variance \(_^\) was estimated. Correlations between variances \(_^\) were estimated in the later stages of model development.

The differences between the observed and individual model-predicted values were modeled as random quantities assumed to follow a normal distribution \(_\sim N\left(0,_^\right)\) for individual \(i\) and the observed value \(j\) of the dependent variable (biomarker) \(k\), where the residual unexplained variability (RUV) \(_^\) was estimated. For each new dependent variable \(k\) added to the model, a flexible combined additive and proportional residual error model was initially used and later challenged for reduction to either an additive or proportional error model.

Parameter estimation method

First-order conditional estimation (FOCE) with interaction was initially used, as it is relatively fast with rich data.Footnote 1 When the model complexity increased, the Stochastic Approximation Expectation–Maximization (SAEM) estimation method was used because it handles complex PKPD models more efficiently [16]. In most models during development, SAEM was preceded by an exploratory Iterative Two-Stage (ITS) estimation step, followed by an Importance Sampling (IMP) evaluation step to evaluate the objective function valueFootnote 2 (OFV). When using SAEM, a full variance–covariance matrix was estimated ($OMEGA BLOCK), and parameters without inter-individual variability were fixed to have a 15% coefficient of variation ($OMEGA 0.0225 FIX), which allowed the SAEM algorithm to search the parameter space efficiently [16]. For SAEM, the AUTO = 1 option with 3000 burn-in iterations (NBURN = 3000) was typically used, followed by 1000 accumulation iterations (NITER = 1000). These settings were tweaked by visually inspecting the convergence of the OFV during the accumulation phase. If the OFV barely changed after a certain number of accumulations, NITER was reduced to a number just above that to save time. Mu-referencing was used whenever possible to achieve greater estimation efficiency [16]. The built-in convergence tester (from the AUTO = 1 option) for SAEM was used to determine when the burn-in phase statistically converged. Additionally, the SEED and RANMETHOD = 3S2P options were used to reduce stochastic noise and thus increase the reproducibility of both SAEM and IMP estimation methods.

The standard errors of the parameter estimates were obtained using the observed Fisher information matrix from the NONMEM covariance step ($COVARIANCE MATRIX = R).

Model qualification criteria

The selection of a candidate model was guided by scientific plausibility as well as numerical and graphical assessments [17].

The OFV was used for comparison between hierarchical models. A nominal p-value lower than 0.01 was chosen to be statistically significant for one degree of freedom—corresponding to a reduction in the OFV of more than 6.64 points. When comparing non-hierarchical models, the Akaike information criteria (AIC) was used [18].

The relative standard errorFootnote 3 (RSE) of an estimated parameter was also used as a model selection criterion. An RSE of < 40% for fixed effects and < 50% for random effects was considered adequate precision. The condition number was also evaluated, with a condition number above 1000 considered indicative of overparameterization.

An observation with a CWRES > 5 was considered an outlier. If an outlier was encountered, sensitivity analysis was performed to evaluate its influence on model parameters.

The ability of the population model to describe the data was assessed using prediction- and simulation-based diagnostics. Various population-based and individual-based goodness-of-fit (GOF) plots were used to evaluate the structural, RUV and IIV-model throughout assessment of model candidates [17]. Visual predictive checks (VPCs) were used to evaluate the predictive performance of the models through stochastic simulations. For all observations, 200 predictions were simulated, resulting in a distribution of predictions. The prediction distribution was then compared with the observations. Because rich data were available, a 95% prediction interval and a 95% confidence interval were used.

Simulations

The developed model was used to stochastically simulate the consequences of an inhibitory drug-target effect of hypothetical compounds upon interaction with specific pathways in the erythropoiesis system. Each drug-target interaction defines a distinct mechanism occurring at different stages of erythrocyte development or hemoglobin synthesis. A treatment period of 120 days was assumed. The doses corresponded to either 0.5 or 2 times the AUC50 of bitopertin, where AUC50 is the area under the concentration–time curve associated with 50% of the maximal pathway inhibition (Imax). The exposure levels translate to approximately 20% and 40% inhibition of the pathway (partial inhibition), respectively [9]. During this process, four potential drug-target effect mechanisms were examined. These mechanisms were:

Mechanism A: Inhibition of hemoglobin synthesis (\(_,\;\text}\)), resembling the effect observed for bitopertin.

Mechanism B: Inhibition of reticulocyte precursor production (\(_,\;\text}\)).

Mechanism C: Inhibition of the differentiation of precursors to reticulocytes (\(_}\)).

Mechanism D: Inhibition of hemoglobin synthesis (\(_,\;\text}\)) and full inhibition of Hbtot-driven feedback.

Additional information on the simulated mechanisms can be found in Supplementary A1.

Software

Model development was performed using NONMEM (v7.4.3) [19], facilitated by PsN (v4.9.3) [20] and Improve (v2.5.1–5); additional simulations were performed in the mrgsolve R package (v1.0.3). For graphical analyses, R Statistical Software (v4.1.2) [21] was used, operating in the RStudio environment (v1.4.1717) [22], on a system running Windows 10 Enterprise.