Quantitative Systems Pharmacology Modeling of Avadomide-Induced Neutropenia Enables Virtual Clinical Dose and Schedule Finding Studies

Avadomide is a cereblon E3 ligase modulator and a potent antitumor and immunomodulatory agent. Avadomide trials are challenged by neutropenia as a major adverse event and a dose-limiting toxicity. Intermittent dosing schedules supported by preclinical data provide a strategy to reduce frequency and severity of neutropenia; however, the identification of optimal dosing schedules remains a clinical challenge. Quantitative systems pharmacology (QSP) modeling offers opportunities for virtual screening of efficacy and toxicity levels produced by alternative dose and schedule regimens, thereby supporting decision-making in translational drug development. We formulated a QSP model to capture the mechanism of avadomide-induced neutropenia, which involves cereblon-mediated degradation of transcription factor Ikaros, resulting in a maturation block of the neutrophil lineage. The neutropenia model was integrated with avadomide-specific pharmacokinetic and pharmacodynamic models to capture dose-dependent effects. Additionally, we generated a disease-specific virtual patient population to represent the variability in patient characteristics and response to treatment observed for a diffuse large B-cell lymphoma trial cohort. Model utility was demonstrated by simulating the avadomide effect in the virtual population for various dosing schedules and determining the incidence of high-grade neutropenia, its duration, and the probability of recovery to low-grade neutropenia. Supplementary Information The online version contains supplementary material available at 10.1208/s12248-021-00623-8.


INTRODUCTION
Neutrophils are a major class of white blood cells (1). Neutrophils mature in the bone marrow, move to and reside in peripheral blood circulation, and migrate to inflamed tissue sites when necessary (2). Here, neutrophils can degranulate, phagocyte microbes, or release cytokines to amplify inflammatory response (3). The blood count of neutrophils (absolute neutrophil count or ANC) is a clinical metric for individual capability to fight infections. Neutropenia is a state of low ANC (4,5), which can occur due to genetic disorders (e.g., cyclic neutropenia) and immune diseases (e.g., Crohn's disease) or may occur as a drug-induced toxicity (6).
IMiDs (immunomodulatory drugs) and CELMoDs (cereblon E3 ligase modulation drugs) are a class of compounds therapeutically active against a number of malignancies. These therapeutics include thalidomide, lenalidomide, pomalidomide (7), and others currently in clinical development (e.g., iberdomide (8)). IMiD/CELMoD compounds bind to cereblon (CRBN) and modulate the affinity of the cereblon E3 ubiquitin ligase complex (CRL4 CRBN ) to its substrates, thereby favoring their recruitment, ubiquitination, and subsequent proteasomal degradation. Avadomide (CC-122) is a novel CELMoD being developed for patients with advanced solid tumors, non-Hodgkin lymphoma (NHL), and multiple myeloma (MM) (9). While research continues towards full elucidation of avadomide activity, it is known that avadomide drives CRL4 CRBN interaction with two hematopoietic zinc finger transcription factors (Ikaros (IKZF1) and Aiolos (IKZF3)) inducing their degradation. These transcription factors are known to promote immune cell maturation (10) and normal B-and T-cell function (11). Avadomide administration is associated with a potent antitumor effect and stimulation of T and NK cells in diffuse large B-cell lymphoma (DLBCL) patients (12).
In a recent phase I trial for avadomide in patients with advanced solid tumors, NHL, or MM (trial identifier: NCT01421524), 85% of patients experienced treatment-emergent grade 3/4 adverse events, primarily neutropenia, followed by infections, anemia, and febrile neutropenia (13). Clinical management of neutropenia includes adjunct therapies to stimulate neutrophil production (e.g., administration of granulocyte-colony stimulating growth factor (G-CSF) as filgrastim), dose reduction, or treatment discontinuation. Another approach to manage avadomide-induced neutropenia is the introduction of an intermittent dosing schedule. For example, 5 days on-followed by 2 days off-treatment (5/7 schedule) improved tolerability and reduced frequency and severity of neutropenia, febrile neutropenia, and infections (13).
In this context, quantitative systems pharmacology (QSP) modeling offers opportunities for in silico exploration of alternative dose and schedules that maximize drug exposure while allowing for toxicity management.
Published models of drug-induced neutropenia (14)(15)(16) are inappropriate for avadomide due to the different mechanisms of action. Avadomide-induced degradation of Ikaros leads to myeloid maturation arrest at a promyelocyte stage and does not affect proliferative neutrophil precursors (13). Hence, we develop a QSP model to represent avadomideinduced neutropenia.
To our knowledge, this is the first model specifically developed for neutropenia caused by block in neutrophil maturation. The model is applied to predict the incidence and the severity of neutropenic events in a virtual DLBCL population across a range of dosing schedules. Such a QSP tool is needed because CELMoDs are a large and growing family of compounds and many CELMoDs developed to date share similar patterns of toxicity.

METHODS
This section details technical and methodological aspects of model development and evaluation. Model development used in vitro neutrophil maturation and clinical ANC data, whereas model evaluation of neutropenia pattern used simulation-generated data.
The model is ordinary differential equation (ODE) based and was integrated using Matlab R2020a ODE routines (32). For model fit we applied the optimization routine fminsearch (33) to minimize an objective function consisting in the weighted sum of absolute normalized difference between model simulation and experimental data.

Model Identifiability and Sensitivity Analysis
Structural identifiability verifies that, given the proposed model structure, it is possible to regress a unique set of model parameters (globally or locally) under the hypothesis of ideal data (noise-free and continuously sampled) (34). This test was conducted in Matlab using the GenSSI 2.0 package (35)(36)(37).
Sensitivity analysis (SA) allows exploration of model input-output structure and supports model development. Global SA (GSA) enables a broad exploration of parameter space. We adopted a Monte Carlo-based method as described in (38) (Supplementary Material 1.1).

Virtual Patient Population
To represent the heterogeneity of ANC data observed in the clinical trial, we generated virtual patient cohorts. A virtual patient consists of a neutrophil life cycle model for which selected parameters are assigned from disease-specific (e.g., glioblastoma (GBM) or DLBCL) probability functions. To obtain parameter empirical distributions, the model is repeatedly fitted to individual clinical ANC data. These distributions are tested for normality by applying the Anderson-Darling test (adtest, Matlab) and smoothed adopting a kernel density estimation (ksdensity, Matlab).

Model Validation
For validation, the model simulations were compared to clinical datasets that were not used during the virtual population development. The comparison was based on a two-sample Kolmogorov-Smirnov (K-S) test. This statistical test determines if the empirical distributions of two sample sets belong to the same distribution. Here, the two sample sets are the model generated ANC and clinical ANC taken at the same time after avadomide administration. This test was executed in Matlab using the kstest2 function.

Quantification/Assessment of Neutropenia Severity
The final goal of the simulation is the quantification of neutropenia incidence for a given avadomide dosing schedule in a virtual patient population. We focused on neutropenia and did not develop an efficacy-pharmacodynamic (PD) model for tumor suppression. We adopted drug exposure (e.g., area-under-thecurve or AUC in the central compartment of the PK model) as surrogate endpoint for efficacy, assuming direct proportionality between exposure and efficacy. Drug exposure was compared to neutropenia severity based on the following clinical parameters: (i) toxicity event (i.e., occurrence of any neutropenic event), (ii) 7-day toxicity event (i.e., neutropenic event lasting for at least 7 consecutive days), (iii) recovery from neutropenia (i.e., recovery to grade 1, meaning at least one ANC measure above grade 2 threshold after a toxicity event), (iv) time to recovery (i.e., time between first toxicity onset and first subsequent ANC above grade 2), and (v) time to ANC nadir. The toxicity events considered were neutropenia grade 3 (ANC below 1E9 neutrophil/liter) and grade 4 (ANC below 5E8 neutrophil/liter). The evaluation of 7-day neutropenia was preferred since grade 4 neutropenia lasting 7 days or more is a dose-limiting toxicity by protocol. Simulation analysis was limited to the first treatment cycle (28 days).

Neutrophil Life Cycle Model Captures Main Stages of Neutrophil Maturation
The QSP workflow is shown in Figure 1. It integrates three modules (i.e., PK, PD, neutrophil life cycle) and accessory operations (e.g., definition of virtual patients, model validation).
The neutrophil life cycle model (Eqs. 1-8) describes the neutrophil formation and maturation processes in bone marrow hematopoietic space, egress to peripheral blood circulation, and terminal death. The model consists in a proliferation pool (proliferation), with proliferation rate k prol ; a sequence of transit stages (transit 1, 2, 3, representing progressive maturation according to in vitro studies (39,40)), with rate constants k tr1 , k tr2 , k tr3 , k tr4 ; a bone marrow reservoir pool (reservoir) of mature neutrophils and final release, with k out rate constant, to peripheral blood (circulation). Circulating neutrophils are subjected to terminal death based on k elim rate, while maturing neutrophils undergo apoptosis based on k d rate constant.
Avadomide reversible and incomplete block of cell maturation in vitro occurs primarily at the late maturation stages of neutrophil development and does not affect proliferative neutrophil precursors (13,(39)(40)(41); therefore, it was applied to the transfer rate between transit 2 and transit 3 (Effect CC-122 in Eqs. 3, 4, and 9). To represent a nonlinear response of maturing neutrophils to drug induce-perturbation of maturation process, k tr,3 expression was modified into a Michaelis-Menten-based functional form (k tr;3 ¼ Vmax KMþTransit2 , see Eqs. 3-4 and Supplementary Materials 2.1 for details). The model includes two regulatory feedback mechanisms of neutrophil maturation: Feedback proliferation (Eq. 7) modulates the proliferation rate based on transit 2 level, and feedback egress (Eq. 8) regulates egress of neutrophils from reservoir pool to peripheral blood. Both feedback mechanisms have a similar functional form, and the exponents (γ and β) modulate the velocity of the control action. For full details of the model formulation refer to Supplementary Materials 2.1.
Avadomide PK and PD Models The avadomide PK is described by a two-compartment PK model (42). The avadomide PD model (Eq. 9) determines the magnitude of neutrophil maturation block as a function of avadomide concentration. PK/PD model details in Supplementary Materials 2.2.

Clinical Trial Data Show High Inter-and Intra-disease Cohort Variability in Longitudinal ANC Patterns
We conducted a preliminary data analysis to explore patterns of longitudinal ANC profiles for the first treatment cycle ( Figure 2) across and within disease cohorts and dosing groups. Because the representation of G-CSF was out of model scope, individual ANC subsequent to the first G-CSF administration were removed prior to further data/model analyses (details in Supplementary Material 2.3). This analysis revealed a significant variability in the longitudinal ANC profiles that associated with both initial patient characteristics (e.g., baseline ANC measures from~2E9 to 8E9 cell/liter, Figure 2a) and treatment dosing schedules (nadir depth normalized to baseline varies within the same disease cohort for different dosing schedules, Figure 2c). These results emphasize the need to generate disease-specific models and the importance of capturing patient variability within individual cohorts.

Model Parameterization Explains Disease Cohort Differences in ANC Patterns
Model parameterization involved a combination of literature information, experimental observations, calculation, and regression.
Because the neutrophil life cycle model (detailed in Supplementary Material 2.1) has a unidirectional and sequential transit compartment structure, most of the parameters can be calculated given one of these transit rates. We informed k elim from literature and fixed k d to a minor/negligible rate (as detailed below), and back-calculated k out , k tr4 , k tr3 , k tr2 , k tr1 , k prol under the assumption of homeostasis (i.e., cell count remains constant in all compartments). Calculation details are shown in Table I.
The half-life of circulating neutrophils in humans is a subject of discussion. Several publications report contrasting data (43)(44)(45)(46), proposing that half-life could range from a few hours to several days. Difficulty in measuring this parameter depends mostly on the cell-labeling system adopted and to the fact that neutrophils can relocate to marginated sites, thereby affecting apparent circulating half-life estimates. Furthermore, neutrophil life-span can change under nonhomeostatic conditions (46). In particular, Dale et al. (47) reported that under neutropenic state, neutrophil life-span doubles (t 1/2 = 9.6 h control vs 20.3 h neutropenia state). Given this knowledge and because the majority of papers report half-life ranging from 4 to 18 h (46), with a recent report measuring 3.8 days (48), we choose a typical value of 15 h, and we double it to 30 h in agreement with enhanced life-span for neutropenia disease state. Finally, because all transit parameters are related, the choice of a different t 1/2 within this range would not lead to significant changes in model outputs.
For initial cell count in the model compartments, because it was not possible to determine neutrophil cell concentration in the human hematopoietic tissue in vivo, we adopted the same approach of Friberg et al. 2002 (24) and fixed the initial cell level in all compartments (excluding the reservoir component) to the initial neutrophil concentration in blood.
The remaining parameters were regressed or fixed to constant values. Regressed parameters include the exponent of the feedback proliferation function (γ); the initial cell level in the reservoir pool (expressed as the ratio of cell level in the reservoir pool divided by cell level in circulation, or Ratio Reserv0/Circ0 ), and K M (in the following expressed as a fraction of the initial cell level in transit 2 compartment, or K M, fraction ). These parameters allow modulation of neutropenia patterns in different disease cohorts (e.g., GBM or DLBCL patients) or across individual patients and are discussed below. Fixed parameters are k d and β. k d was introduced above as a maturing cell death rate. The in vitro maturation assay showed that avadomide induces a reversible maturation block with no significant change in cell viability. However, apoptosis of maturing cells is a biologically recognized process, and it is possible to speculate that in vivo neutrophils undergoing long-term maturation block may experience enhanced apoptosis. Based on this, we The neutrophil life cycle model is based on a compartmental structure. The proliferation pool represents committed proliferative neutrophil precursors. From a model idealization standpoint, these cells have specific characteristics: they can proliferate but not self-renew and can proceed to subsequent maturation stages, represented in the model as a sequence of transit compartments. These compartments (i.e., transit 1, transit 2, and transit 3) do not have a direct biological counterpart but here are intended to capture the fact that progressive maturation implies a time delay, in line with previously published implementations of neutrophil maturation models. Once maturation is completed, cells are stored in a bone marrow reservoir pool, awaiting egress into peripheral blood circulation. Circulation pool represents circulating neutrophils (i.e., level of neutrophils in blood, comparable to clinical ANC). Finally, circulating neutrophils are subjected to terminal elimination (cell death) included this process in the model with an arbitrarily assigned small rate (i.e., 0.001 h -1 or~4% of k tr maturation rates departing from the same compartments). The parameter β controls egress rate from the bone marrow reservoir pool. The biological mechanism controlling neutrophil egress from bone marrow is complex and only partially understood (49). We fixed β to a high value based on the clinical observation that, even in the presence of avadomide block, circulating ANC was maintained at a baseline level for several days despite compromised bone marrow maturation, suggesting that the egress of mature neutrophils from bone marrow is sustained and prompt.
The model was initially fitted to data from GBM patients. Those patients did not receive previous lines of bone marrow depleting treatments and therefore represent the closest match to a healthy bone marrow condition before avadomide treatment. The model was fitted simultaneously to all GBM dose groups in order to regress a single parameter set representative of the GBM patient population (Figure 3a). At this step, five parameters were fitted. Two of those parameters are PD specific (EC 50,PD and n PD ) and three are disease group specific (γ, Ratio Reserv0/Circ0 , K M, fraction ). Once regressed, PD parameters were kept constant for any other avadomide simulation/fit under the assumption that drug effect is reproducible across the disease cohorts. The three disease group-specific parameters were instead refitted per disease group, because these parameters are representative for the bone marrow state and thus change across disease cohorts.
For model fit to the DLBCL median profiles (i.e., graydotted lines in Figure 3b), the parameters γ, Ratio Reserv0/Circ0 , and K M, fraction were regressed using the GBM parameter values as an initial guess. This operation served multiple purposes: (i) determine typical parameter values of DLBCL patients, (ii) explore whether parameter value differences between GBM and DLBCL could explain biological differences between the two patient groups, and (iii) determine initial parameter estimates for the subsequent step of patient-specific model fits. Figure 3b shows a model fit to median DLBCL ANC data, and Table I   (i.e., Ratio Reserv0/Circ0 ), the extent of proliferative response to avadomide maturation block (i.e., γ), and idiosyncratic capacity to contrast maturation block (i.e., K M,fraction ) are reduced in DLBCL compared to GBM.

Virtual Patient Cohort
Four model parameters allow for characterization of individual patients: ANC level at baseline, Ratio Reserv0/Circ0 , K M,fraction , and γ. Briefly, the ANC level at baseline is the neutrophil count in blood before treatment start. Ratio Reserv0/ Circ0 is the individual initial level of mature neutrophils stored in the bone marrow. K M,fraction regulates changes to neutrophil transfer from transit 2 to transit 3 when transit 2 cell level deviates from its homeostatic value. γ controls the magnitude of proliferative response to the avadomide-induced perturbation of neutrophil maturation.
Starting from the DLBCL reference parameter set, the model was refitted to individual ANC profiles in the DLBCL cohort, thereby generating a set of values for each parameter. Because not all parameter value distributions are normal, we kept the parameter empirical distributions as they are (i.e., without replacing them with parametric models) and adopted kernel density estimation to estimate the probability density function (Figure 4a).
Finally, virtual patients were created by independent random sampling from the parameter value probability distribution functions (parameter values are assumed independent, meaning that there is no conditional probability for parameter values given the value of other parameters). The virtual cohorts generated for this analysis included 1,000 virtual patients (Figure 4b).

Model Identifiability and Global Sensitivity Analyses
The model was tested for identifiability considering the three individualized parameters (γ, K M,fraction , Ratio Reserv0/ Circ0 ) and specifying that observations are only available for the circulation compartment. K M,fraction , and Ratio Reserv0/Circ0 are globally structurally identifiable, while γ is locally identifiable.
We used GSA to rank parameters by importance in determining changes to the simulated ANC profile (full results in Supplementary Materials 2.4). GSA results support the choice of γ and Ratio Reserv0/Circ0 as individual parameters for the generation of the virtual patient population and indicate that K M,fraction is likely to contribute poorly toward differentiating virtual patients. For the present application, we acknowledge the minor role of this parameter, which could nonetheless be relevant for model application in the context of other indications, and it is therefore kept in the virtual patient generation workflow.

Virtual Population of DLBCL Patients Reproduces Clinically Observed Longitudinal ANC Profiles
The virtual DLBCL patient population was validated by simulating the same treatment received by two clinical trial cohorts (avadomide 3 mg on a 5/7 and QD schedule, data not used to generate the virtual population) and then testing  Figure 5 shows how these distributions were found being equivalent at all tested times for the 3-mg QD group and for 4 of 5 times for the 3-mg 5/7 group. Qualitatively, simulated longitudinal ANC for the virtual cohort ( Figures 5, 6, and 8) compare nicely to clinical data ( Figure 3): ANC levels are quite stable until day 8-12 and then drop until about day 20, followed by stable low ANC count with administration schedule-driven fluctuations.

Model is Applied to Explore Doses and Schedules
Avadomide administration to the virtual DLBCL cohort (1000 virtual patients) was simulated for all combinations of 7 doses (i.e., 2, 3, 4, 5, 6, 7, 8 mg) and 6 schedules (i.e., 3/7, 5/7, 7/14, 14/28, 21/28, 28/28), totaling 42,000 simulations. Next, individual predictions of ANC profiles were processed to determine whether or not avadomide caused grade 3 or 4 neutropenia, its duration, the recovery, and the time to recover. Collective analysis determined the percentage of patients expected to experience toxicity and possibly recover from it within the first drug administration cycle. Here we report a selection of representative results (full results in Supplementary Materials 2.5). Figure 6 shows the longitudinal ANC profiles for the same virtual cohort receiving 6 mg of avadomide on the 5/7 or 21/28 schedule. In terms of exposure, the two schedules allow similar total dosing and PK exposure over the first cycle (20 doses and 1417 ng/ml*h AUC cycle1 vs 21 doses and 1515 ng/ml*h AUC cycle1 , for schedules 5/7 and 21/28, respectively). Simulations show that until exhaustion of the reservoir pool, the ANC level remains stable, whereas at later time points (typically after day 10 post administration) ANC start dropping towards neutropenic levels. The schedule 5/7 shows   Table II shows the incidence of high-grade neutropenia and recovery for (i) different schedules at the same dose (4 mg) and for (ii) the same schedule at different doses (5/7, 2 to 8 mg).
Based on Table IIA, drug exposure (measured as AUC) increases with the total number of dosing days while C max Figure 5. Model validation results. a Avadomide 3-mg QD. Top: longitudinal ANC profiles, virtual cohort (1000 subjects) = gray-solid, clinical cohort (18 patients) = blue-dotted. Bottom: K-S test for equivalence of cumulative distribution profiles (with 5% significance level P value ). b Avadomide 3mg 5/7 day. Top: longitudinal ANC profiles, virtual cohort (1000 subjects) = gray-solid, clinical cohort (14 patients) = blue-dotted. Bottom: K-S test for equivalence of cumulative distribution profiles (with 5% significance level P value ). Virtual and clinical ANC distributions were taken at day 1, 8, 16, 22, and 28 and compared using the two-sample K-S test. Distribution equivalence rejected only for 3mg 5/7 at day 22 (i.e., equivalence verified at day 1, 8, 16, 28, but not at day 22). The K-S test confirms that all the simulated profiles for the virtual cohort run to a clinically plausible state increases with the number of consecutive dosing days. For neutropenia, the incidences of both grade 3 and 4 neutropenic events increase with consecutive dosing days, with the exceptions of 5/7 which shows a slightly higher incidence than 7/14. In contrast the incidence is not directly dependent to the total dose received, as shown by the differences between 7/14 vs 14/28 or 5/7 vs 21/28. Interestingly, the incidence of grade 3 and 4 events is very similar for schedules 21/28 and 28/28. However, this similarity is not found for neutropenia maintained for at least 7 consecutive (7+) days, where we observe a substantial difference between schedules 21/28 and 28/28, which show the incidence of 36.6% and 45.6% (for grade 3, 7+ days), respectively. For 28/28 single and 7+ day, neutropenia has the same total incidence, while intermitted schedules show a reduction of 7+ neutropenic events compared to single events. In terms of recovery, all the intermittent schedules with at least 7 days of dose interruption show substantial recovery (i.e., 66% (12.5/19), 83% (28/ 33.7), and 84% (38.5/45.4) of virtual patients that experienced neutropenia grade 3 recovered above grade 2 for 7/14, 14/28, and 21/28, respectively). In contrast, no recovery was determined for 3/7 and 5/7 schedules. For schedules that allow recovery, the recovery time increases nonlinearly with consecutive dosing days (i.e., 4.7, 6.3, and 11.2 days were necessary on average to recover from grade 3 to above grade 2 for schedules 7/14, 14/28, and 21/28, respectively). Table IIB focuses on the 5/7 schedule: both AUC and C max increase linearly with the dose, the incidences of both grade 3 and 4 neutropenic events increase with dose, and recovery is absent or minimal at all doses. Figure 7 shows a bar plot comparison of toxicity and recovery across schedules for two doses (4 or 6 mg). Bars are schedule-specific and are ordered by increasing drug exposure. The percentage of patients experiencing toxicity increases with the number of consecutive dosing days. This pattern is not verified for 5/7 vs 7/14 likely because of the combined effect of a similar number of dosing days (5 vs 7 days) and the difference in the dosing holiday duration (2 vs 7 days). Recovery from grade 3 is substantial (>80%) and very similar for 14/28 and 21/28 and increases with dose for schedules 7/14 and 14/28, but not for 21/28. Increase in dose from 4 to 6 mg associates with higher recovery from grade 4. Schedule 5/7 shows some lower toxicity compared to other schedules but offers little or no recovery. Figure 8 shows the time of nadir for five different schedules. Schedule 5/7 shows a bimodal time of nadir with 9% of patients having nadir at day 20 and~91% at day 27. Schedule 7/14 and 21/28 show nadir at day 21, consistently with the start of the latest dosing holiday for cycle 1. Schedule 14/28 shows nadir in the interval of day 15 to 17. Finally, daily dosing (schedule 28/28) results in a progressive increase of the virtual patients having ANC nadir in the interval of day 21 to day 28.

DISCUSSION
In this paper, we have presented a QSP model for avadomide-induced neutropenia. We applied this model to virtually explore the pattern and the incidence of neutropenia across dosing schedule scenarios in a DLBCL patient population treated with avadomide. Model development followed good practice standards as described in Bai et al. 2019 (17).
The neutrophil life cycle model developed describes neutrophil maturation and transit stages from bone marrow to peripheral blood and captures the avadomide-specific mechanism of induction of neutropenia. Since this mechanism is different from chemotherapy-induced neutropenia, published models (such as the Friberg model (24)) could not be applied to address the needs of our study. A major difference of our model compared to the Friberg model (24) is that proliferation rate is not controlled by ANC-level changes compared to baseline in peripheral blood. The latter mechanistic implementation was not well-suited for the description of the CELMoD-driven neutrophil maturation block and caused indefinite accumulation of neutrophils at the maturation blocked stage and excessive proliferation (because during maturation block, proliferation was continuously stimulated by the sub-baseline ANC level). Additionally, a first-order modeling of the cell transit through maturation stages is not suitable for CELMoD-like maturation The ANC baseline distribution (i.e., ANC at t=0) is the same because the same virtual patients are simulated for both dosing schedules. The two schedules enable very similar PK exposure over the first treatment cycle; however, the neutropenia pattern is quite different: schedule 21/28 shows deeper ANC drop and protracted toxicity, followed by strong recovery once the treatment is interrupted. In contrast, schedule 5/7 offers a mitigated incidence of high-grade toxicity, with only limited recovery during dose interruption  Figure 7. Bar plot analysis for toxicity and recovery for different schedules at 4 mg (a) and 6 mg (b). Grades 3 and 4 single indicate the percentage of virtual patients experiencing at least one event of neutrophil level below the respective toxic threshold. Grades 3 and 4 7 days indicate percentage of virtual patients experiencing an extended and uninterrupted toxicity for at least 7 days. Recovery Gr3 to above Gr2 and Gr4 to above Gr2 indicate the percentage of patients that recovered to grade 1 (i.e., above grade 2) relative to the patients that experienced toxicity. This analysis is limited to the first treatment cycle block. For example, the first-order-based transit (i.e., transit rate constant*cell level in the upstream compartment) in presence of CELMoD-effect (i.e., reduction of transit rate constant) causes an increase of cell level, which over time would compensate for transit rate constant reduction and ultimately mask drug effect. Accordingly, we adopted a Michaelis-Mentenlike function for transit stage 2 which allowed an asymptotic behavior of the flow out of transit 2 despite an increase in accumulated maturing neutrophils.
In terms of the workflow, the clinically observed variability of ANC supported extending model simulation from a single median virtual patient to a virtual patient population. The DLBCL virtual cohort utilized in our simulations was validated comparing the cumulated distributions of the clinical and the virtual cohorts ANC at selected time points. This approach allowed for both qualitative and quantitative evaluation of equivalence of the two empirical cumulated distributions. An alternative and commonly adopted approach, like the visual predictive check, is conceptually similar in terms of comparing virtual vs clinical distributions, but it is more qualitative in nature.
The heterogeneity of the virtual population is observable in the simulated ANC profiles in terms of initial baseline, neutrophil reservoir pool size (ANC starts dropping from baseline level at different times), and idiosyncratic variability in response to maturation block (visible as an overlapping profile in the recovery time interval). A limitation of the current implementation is that population PK was not included, as that would improve the representation of the variability across the virtual population.
Model utility was demonstrated by simulating avadomide administration to a virtual DLBCL cohort. Since it was not possible to develop an avadomide efficacy module in the absence of specific biomarkers or tumor suppression data, the drug exposure (i.e., AUC in central PK model compartment) was used as a reference to contrast schedule toxicity.
Simulation results address different aspects of neutropenia pattern modulation by choice of dosing schedule. Frequent dosing (i.e., schedules 28/28 and 5/7) produces high systemic exposure along with the highest incidence of neutropenia. The 2-day dosing holiday on the 5/7 schedule is sufficient to reduce significantly the total incidence of neutropenia in the virtual population (e.g., at the 4-mg dose, the schedule 5/7 compared to 28/28 gives~28% less exposure, but it lowers the incidence of neutropenia grade 3 by~44%). However, a 2-day holiday does not grant measurable recovery from high-grade neutropenia. This suggests that for avadomide in DLBCL patients a longer dosing holiday should be considered in case a more substantial recovery is desired. For example, compared to 5/7 and 28/28, all other tested schedules with a measurable incidence of neutropenia enable substantial recovery (Figure 7). It is noted that the exploration of neutrophil recovery rate during dosing holiday is only possible with model-based tools since trial patients are typically undergoing sequential cycles of treatment and receive concomitant medications for the mitigation of neutropenia.
Regarding the analysis of high-grade neutropenia lasting at least 7 consecutive days (7+ day), among those schedules allowing dosing interruption (excluding 28/28), schedule 21/28 results in a higher incidence of prolonged neutropenia, coherently with the 21-day continuous dosing not allowing for intermittent recovery. The schedule 5/7, despite some mitigation enabled by the 2 days of dosing interruption, produces a 7+ day neutropenia comparable to schedule 14/28. Schedule 7/14 shows the best performance in terms of minimizing 7+ day toxicity at dose level 4 to 6 mg. Further, results show that under continued dosing, the maximal neutropenia would be reached by day 21 (or a few days earlier), since the total incidence of high-grade neutropenia is nearly equivalent for schedule 21/28 and 28/28 (Table II). Finally, the model enables predictions of the time at which the most severe neutropenia is reached (i.e., ANC nadir, Figure 8), showing that nadir time is primarily controlled by the schedule of choice, rather than the dose level.
Collectively, these model-based results show that the choice of dose and schedule offers a powerful handle to modulate the neutropenia incidence, duration, and recovery in the patient population. These results demonstrate the model potential applicability as a support tool to inform decision-making in the clinic, by informing the dose and schedule that would maximize drug exposure while allowing for toxicity management. The calibration of this model to a different compound (with similar MoA) and/or to a different indication would allow the in silico exploration of doses and schedules. Simulation results should be interpreted in the light of clinical protocol definitions for doselimiting toxicity and maximum tolerated dose as well as efficacy considerations.
We aimed at providing a detailed explanation of all the assumptions that went into the development of the model, in order to make it easier for modelers and scientists in general not only to reproduce it using their favorite coding language (50) (R code available in the Supplementary Materials) but also to expand and leverage it for the development of models of neutropenia that are applied to other diseases (51,52) and therapeutic contexts (e.g., beta-hemoglobinopathies (53)).

CONCLUSIONS
Neutropenia is a major treatment-emergent and doselimiting toxicity in trial patients treated with avadomide. Intermittent dosing is an option to manage this toxicity and different combinations of dose and schedule enable controlling the toxicity-efficacy tradeoff. Here we presented a QSP model for avadomide-induced neutropenia, which includes a mechanistic model of neutrophil life cycle combined with avadomide PK and PD. The complete workflow allowed to capture the disease cohort variability and enabled performing simulations for several dosing schedule scenarios, aiming at screening options that would minimize neutropenia while enhancing drug exposure.
This model is the first one developed specifically for neutropenia caused by block in neutrophil maturation, informed by preclinical and clinical data (13,(39)(40)(41). We anticipate further opportunities to apply and develop and demonstrate the relevance of this model given the potential use of avadomide and other CELMoD compounds, either as single agents or in combination to treat a range of indications.