A systems pharmacology model for gene therapy in sickle cell disease.

We developed a mathematical model for autologous stem cell therapy to cure sickle cell disease (SCD). Experimental therapies using this approach seek to engraft stem cells containing a curative gene. These stem cells are expected to produce a lifelong supply of red blood cells (RBCs) containing an anti-sickling hemoglobin. This complex, multistep treatment is expensive, and there is limited patient data available from early clinical trials. Our objective was to quantify the impact of treatment parameters, such as initial stem cell dose, efficiency of lentiviral transduction, and degree of bone marrow preconditioning on engraftment efficiency, peripheral RBC numbers, and anti-sickling hemoglobin levels over time. We used ordinary differential equations to model RBC production from progenitor cells in the bone marrow, and hemoglobin assembly from its constituent globin monomers. The model recapitulates observed RBC and hemoglobin levels in healthy and SCD phenotypes. Treatment simulations predict dynamics of stem cell engraftment and RBC containing the therapeutic gene product. Post-treatment dynamics show an early phase of reconstitution due to short lived stem cells, followed by a sustained RBC production from stable engraftment of long-term stem cells. This biphasic behavior was previously reported in the literature. Sensitivity analysis of the model quantified relationships between treatment parameters and efficacy. The initial dose of transduced stem cells, and the intensity of myeloablative bone marrow preconditioning are predicted to most positively impact long-term outcomes. The quantitative systems pharmacology approach used here demonstrates the value of model-assisted therapeutic design for gene therapies in SCD.

WHAT IS THE CURRENT KNOWLEDGE ON THE TOPIC?

Gene therapy presents a potentially curative treatment for sickle cell disease (SCD). Due to high cost and long timelines of treatment stabilization, it is difficult to predict treatment efficacy and long‐term outcome at early time points post cell transplantation, and challenging to optimize treatment design.

WHAT QUESTION DID THIS STUDY ADDRESS?

We developed a quantitative systems pharmacology (QSP) model to understand how varying specific treatment parameters affects short and long‐term measures of treatment efficacy.

WHAT DOES THIS STUDY ADD TO OUR KNOWLEDGE?

The QSP model provides a rational way of predicting the effect of varying treatment parameters on patient outcomes.

HOW MIGHT THIS CHANGE DRUG DISCOVERY DEVELOPMENT, AND/OR THERAPEUTICS?

Our QSP model and model analysis can guide rational therapeutic design of gene therapy for SCD and other genetic disorders.

INTRODUCTION

Sickle cell disease (SCD) is a genetically inherited blood disorder caused by an abnormal form of hemoglobin, HbS (α2βs 2) due to a single point mutation (E6 V) in the first exon of the β‐globin gene on chromosome 11. Deoxygenated HbS tends to polymerize into long rod‐like structures conferring red blood cells (RBCs) their characteristic sickle‐shaped morphology, reduced oxygen‐carrying capacity, and a short lifespan. These misshapen RBC obstruct capillaries, setting off recurring, painful, vaso‐occlusive crises in patients. Patients with SCD suffer from a host of other symptoms, including hemolytic anemia, organ damage, bacterial infections, and have a lower‐than‐average life expectancy. SCD is largely incurable at this time, other than via allogeneic stem cell transplantation, which is limited by the paucity of suitably matched donors. Current standard treatments primarily seek to alleviate symptoms of the disease.

Several experimental autologous stem cell transplantation therapies for SCD are currently under development. These therapies harvest stem cells from a patient, treat them ex vivo with a self‐inactivating lentiviral vector containing a treatment gene, and transplant the treated stem cells back into the patient. The gene product is either a corrected β‐globin that forms HbA (α2β2), or a γ‐globin that forms fetal hemoglobin, HbF (α2γ2). Preclinical studies show that expression of the γ‐globin gene leads to efficient formation of HbF (α2γ2) in adult erythroid progenitor cells. HbF inhibits sickling due to HbS, , , and HbF induction corrects hematologic and pathologic defects associated with SCD. Viable engraftment of treated stem cells is therefore expected to cure SCD by providing a lifelong supply of RBCs containing anti‐sickling hemoglobin after a single treatment.

Widespread clinical use of gene therapy for SCD has been limited due to the complexity of this multistep procedure. Treatment outcomes, such as RBC counts, and blood hemoglobin levels do not stabilize until several months to more than a year after the initial myeloablation and stem cell infusion. Therefore, it is difficult to predict treatment efficacy and long‐term outcome at the outset, and challenging to optimize treatment design. Motivated by these challenges, we developed a quantitative systems pharmacology (QSP) model of erythropoiesis and autologous stem cell transplantation. Our goal was to understand how varying specific treatment parameters, such as initial dose and ex vivo transduction efficiency, affects short and long‐term measures of treatment efficacy.

Hematopoiesis is the complex process in which many mature blood cells are derived from a single hematopoietic stem cell (HSC), through an ordered hierarchy of intermediate differentiation stages commonly defined by specific combinations of cell surface markers. Previously published mathematical models of hematopoiesis have used stochastic simulations to capture the propagation of a mutant lineage, and to model the successive flux amplification needed for stable hematopoiesis. A published model on gene therapy in SCD used a discrete time age‐structured model to predict the relationship between degree of bone marrow chimerism and the fraction of non‐sickled RBCs. These studies do not explicitly consider hemoglobin dynamics, or the physiology of RBCs in peripheral blood, which are important clinical determinants of treatment success.

Our approach was to construct a multiscale, mathematical model of erythropoiesis and hemoglobin production using ordinary differential equations (ODEs). Alternative parameterizations of this model produce steady states that match healthy and SCD phenotypes on key clinical measures. The model was used to simulate transplantation of autologous stem cells containing an exogenous γ‐globin (γNew) that associates with endogenous α to produce a novel anti‐sickling hemoglobin (HbNew). Treatment simulations predict the dynamics of RBCs derived from transduced stem cells, and the levels of HbNew in peripheral blood. Post‐treatment dynamics in the short and long‐term are shown to be driven by differences in regenerative potential of short and long‐lived stem cells as previously reported. Using a sensitivity analysis, we quantified the effect of varying treatment parameters on various measures of treatment efficacy. This detailed, quantitative analysis is intended to guide rational therapeutic design of gene therapy for SCD and other genetic disorders.

METHODS

Erythropoiesis model

Erythropoiesis is initiated by a self‐renewing LT‐HSC compartment that follows a logistic population growth model. The steady‐state population approaches a set carrying capacity, N. The following ODE model endogenous LT‐HSC, e(t), and transduced LT‐HSC, y(t): for the healthy and untreated SCD phenotypes y(t) = 0 as there are no transduced cells. Here, and are the LT‐HSC population growth and death rate constants, respectively, and is the first order rate constant for differentiation of LT‐HSC into ST‐HSC. These are computed using: where τ LT is the mean residence time in the LT‐HSC compartment.

Cellular flux between successive differentiation stages is represented by unidirectional arrows in Figure 1a. These are implemented as first order processes. Dynamics of each progenitor compartment downstream of LT‐HSC, namely the ST‐HSC, MPP, CMP, burst‐forming unit‐erythroid (BFU‐E), and CFU‐E, is controlled by two parameters: a fold‐amplification (a) and a mean residence time (τ). The resulting system of ODE and model parameterization for healthy and SCD phenotypes is described in the Supporting Information and Table 1.

FIGURE 1

Multiscale model of erythropoiesis and hemoglobin assembly. (a) Schematic of erythropoiesis model implemented with connected compartments to represent successive differentiation stages from LT‐HSC to RBC. The LT‐HSC compartment is self‐renewing and follows a logistic growth model. ST‐HSC through CFU‐E compartments each amplify the influx of cells by the fold‐change () indicated on top. The mean residence time (τ) in each compartment is indicated below the compartment. The variable amplification in the CFU‐E compartment is driven by the feedback from blood oxygen (c). There is no further cellular replication beyond the CFU‐E stage. Sickled RBC lifespan is assumed to be 10‐fold shorter than a healthy RBC lifespan. (b) Model reactions for hemoglobin assembly that are embedded within the RET and RBC pools. All monomers are synthesized only within reticulocytes. Monomers associate to form heterodimers that further assemble into functional hemoglobin tetramers. Monomers and dimers degrade with first order kinetics. (c) Feedback regulation of RBC production via blood oxygen. Venous blood oxygen level is linked to the variable CFU‐E amplification via an exponential function. The function was parameterized to match known reference values for healthy and SCD phenotypes, indicated by the colored dots. BFU‐E, burst‐forming unit‐erythroid; CFU‐E, colony‐forming unit erythroid; CMP, common myeloid progenitor; LT‐HSC, long‐term hematopoietic stem cell; RBC, red blood cell; RET, reticulocyte; SCD, sickle cell disease; ST‐HSC, short‐term hematopoietic stem cell

TABLE 1

Mean residence time and fold‐amplification within each compartment in the erythropoiesis model

Cell type	Compartment name	Mean residence time, τ (days)	Amplification, a [# of divisions = log 2(a)]	Source
Long‐term hematopoietic stem cell	LT‐HSC	100	NA	44
Short‐term hematopoietic stem cell	ST‐HSC	20	1000 [10]	44
Multipotent progenitor cell	MPP	2	1000 [10]	44
Common myeloid progenitor	CMP	4	16 [4]	Estimated by flux matching (See SI for details)
Burst forming unit‐erythroid	BFU‐E	7	32 [5]	20
Colony forming unit‐erythroid	CFU‐E	7	32 [5]
Reticulocyte	RET	3	NA	45
Red blood cell (erythrocyte)	RBC	120 (healthy) 12 (SCD)	NA	29

Abbreviations: NA, not applicable; SCD, sickle cell disease.

Hemoglobin model

Hemoglobin assembly is modeled using an embedded set of molecular reactions within the RET and RBC pools (Figures 1b and 2b). Globin monomers are synthesized within the RET pool, and these assemble into dimers and tetramers. Hemoglobin assembly reactions continue within the RBC compartment. The Supporting Information lists rate laws for hemoglobin assembly and Table 2 lists the parameters used in nominal simulations.

FIGURE 2

Model of autologous stem cell transplantation gene therapy for SCD. (a) Schematic of the initial transfusion of ex vivo CD34+ cells containing transduced (colored) and non‐transduced (white) cells. The fraction of CD34+ cells of each progenitor type is indicated on the corresponding arrow originating from the initial dose. Fold‐amplification within each progenitor compartment is the same as in the baseline model (Figure 1). Transduced cells are assumed to contain a single copy of the vector for the novel exogenous γ‐globin. The lifespan of transduced cells is a model parameter (see Table 3) (b) New hemoglobin assembly reactions in RET and RBC derived from transduced progenitors. The gene product γNew is produced within RET and dimerizes with endogenous α. αγNew dimers assemble into the novel anti‐sickling hemoglobin HbNew. (c) Simulated preconditioning to deplete endogenous progenitors on day 0. Each progenitor pool from ST‐HSC through CMP is depleted by 90% of the pretreatment baseline. (d) Post‐treatment time course of RBC (top panel, 1B = 109 and 1 T = 1012), hemoglobin levels (middle panel), and venous blood oxygen (bottom panel). (e) Post‐treatment time course of transduced LT‐HSC, ST‐HSC, MPP and CMP (top panel, 10 K = 104 and 1 M = 106), RBC derived from transduced LT‐HSC and from transduced ST‐HSC, MPP, and CMP pools combined (middle panel, orange and green curves, respectively. 1 T = 1012), and the ratio of transduced to endogenous RBC (bottom panel) . BFU‐E, burst‐forming unit‐erythroid; CFU‐E, colony‐forming unit erythroid; CMP, common myeloid progenitor; LT‐HSC, long‐term hematopoietic stem cell; MPP, multipotent progenitor cell; RBC, red blood cell; RET, reticulocyte; SCD, sickle cell disease; ST‐HSC, short‐term hematopoietic stem cell

TABLE 2

Hemoglobin assembly model parameters

Parameter	Symbol	Units	Value	Source
Reticulocyte volume	V_RET	pL (= 10−12 L)	0.09	46
Erythrocyte volume	V_RBC	pL (= 10−12 L)	0.09	47
Blood volume	V_blood	L	5	48
Hemoglobin molecular weight	MW_Hb	g/mol	64,500	49
α globin synthesis rate	k_synα	nmol/day/cell	1.5 × 10−6	Estimated (see SI for details and Figure S1 for sensitivity analysis)
β globin / α globin synthesis rate ratio	ratio_synβ	Unitless	0.5	Assumed based on allele numbers (see SI for details and Figure S1 for sensitivity analysis)
γ globin / α globin synthesis rate ratio	ratio_synγ	Unitless	0.03	49
δ globin / α globin synthesis rate ratio	ratio_synδ	Unitless	0.04	49
Free monomer half life	thalf_monomer	Days	0.25	Assumed (see Figure S1 for sensitivity analysis)
Bimolecular binding on rate constant	k_onBimolecular	1/(nmol.s)	10−5	50
αβ dimer dissociation constant	K_dαβ	nM	10−2 (for αβ^s) 10−3 (for αβ)	(see Figure S1 for sensitivity analysis) 51 , 52
αγ dimer dissociation constant	K_dαγ	nM	10−5
αδ dimer dissociation constant	K_dαδ	nM	10−2
HbS (α_2β2s) tetramer dissociation constant	K_dHbS	nM	100
HbA (α_2β_2) tetramer dissociation constant	K_dHbA	nM	100
HbF (α_2γ_2) tetramer dissociation constant	K_dHbF	nM	100
HbA2 (α_2δ_2) tetramer dissociation constant	K_dHbA2	nM	100

TABLE 3

Nominal treatment parameters

Parameter	Symbol	Units	Value	References
Initial cell dose	Dose	Count	108	41, 42
Mean vector copy number	VCN	Unitless	0.5	Based on Scala et al (2018) see Figure 3 for sensitivity analysis
Fractional depletion of LT‐HSC, ST‐HSC, MPP, and CMP by preconditioning	Depletion	Unitless	0.9	Assumed. See Figure 3 for sensitivity analysis
Fraction of CD34+ cells that are LT‐HSC	f_CD34LT	Unitless	10−6	See 53 Figure 3 for sensitivity analysis
Fraction of CD34+ cells that are ST‐HSC	f_CD34ST	Unitless	10−5	Assumed. See Figure S3 for sensitivity analysis
Fraction of CD34+ cells that are MPP	f_CD34MPP	unitless	10−3	Assumed. See Figure S3 for sensitivity analysis
Fraction of CD34+ cells that are CMP	f_CD34CMP	Unitless	10−2	Assumed. See Figure S3 for sensitivity analysis
Fraction of γNew synthesis relative to total (βs + γNew) globin synthesis	f_synγNew	Unitless	0.5	Assuming same synthesis rates for βs and γNew
αγNew dimer dissociation constant	K_dαγNew	nM	10−5	Assumed to be same as endogenous αγ
HbNew (α2γNew 2) tetramer dissociation constant	K_dHbNew	nM	100	Assumed same as endogenous HbF
Mean lifespan of RBC in transduced branches	τ_RBCtrans	Days	90	Assumed. See Figure 3 for sensitivity analysis

Abbreviations: CMP, common myeloid progenitor; LT‐HSC, long‐term hematopoietic stem cell; MPP, multipotent progenitor cell; RBC, red blood cell; RET, reticulocyte; ST‐HSC, short‐term hematopoietic stem cell.

10−2 (for )

10−3 (for )

Feedback on CFU‐E amplification

The number of divisions in the CFU‐E stage is variable, and driven by blood oxygenation. To implement this feedback, we computed the venous blood oxygen level from the blood hemoglobin levels in the model using: vO2 = (0.74 × [HbA]g/dl + 0.68 × [HbS]g/dl + 0.88 × [HbF]g/dl + 0.88 × [HbNew]g/dl) × 1.34 ml/dl. The coefficient for each hemoglobin is its reported saturation at the partial pressure of oxygen in venous blood, and 1.34 ml/dl is the maximum oxygen carrying capacity of hemoglobin.

The computed venous oxygen level was empirically linked to the CFU‐E amplification using an exponential function: (Figure 1c), motivated by the reported exponential relationship between blood epo level and hemoglobin level. , The exponential was parameterized by anchoring it to the reference values for healthy and SCD phenotypes.

Treatment simulations

The model structure was updated as shown in Figure 2a,b to simulate gene therapy with ex vivo transduced CD34+ cells containing a novel gamma globin (γNew). Only a fraction of cells in the initial dose (colored in blue) carry the treatment vector. Transduced precursor cells populate four new branches in the erythropoiesis model, each originating from a specific precursor cell type, seeded in relative proportions indicated in Figure 2a and Table 3. Non‐transduced cells of each precursor type are assumed to be distributed in the same relative proportions and added to the corresponding endogenous pools. The synthesis of γNew and production of HbNew were implemented within the RET and RBC pools of each transduced branch (Figure 2b and Supporting Information).

Treatment simulations were initiated by first depleting endogenous precursor compartments from LT‐HSC through CMP to 10% of their pretreatment baseline values to simulate preconditioning (Figure 2c). All endogenous cell numbers downstream of CMP, and protein levels were left unchanged from pretreatment baselines. In essence, treatment simulations augment the original model with a new set of rate laws, and perturb the pretreatment baseline by applying preconditioning (depletion of precursor cells) and transplantation (addition of precursor cells).

Sensitivity analysis

For the long‐term analysis, each treatment parameter was varied over a range holding all other parameters fixed at their nominal values. Treatment was simulated for each such parameterization and simulation outcomes were recorded at 2 years post‐treatment (Figure 3). Sensitivity analyses to endogenous model parameters and to treatment parameters that affect short‐term outcomes are described in the Supporting Information.

FIGURE 3

One parameter at a time sensitivity analysis of treatment simulations. Each model parameter listed in the column title is varied over a range of values shown on the x‐axis, with all other parameters held fixed at their nominal values indicated by dots. Predicted outcomes shown on the y‐axis are at 2 years post‐treatment. LT‐HSC, long‐term hematopoietic stem cell; RBC, red blood cell; ST‐HSC, short‐term hematopoietic stem cell; VCN, vector copy number

RESULTS

A multiscale model of erythropoiesis and hemoglobin production

In this study, we developed a mathematical model of erythropoiesis and hemoglobin production based on known biology. , , , Erythropoiesis is initiated by long‐term hematopoietic stem cells (LT‐HSCs) in the bone marrow. These undergo several rounds of differentiation, finally producing reticulocytes (RETs) that are released into the blood where they mature into RBCs. We modeled this process using connected compartments with unidirectional flow to represent the successive differentiation stages (Figure 1a).

Self‐renewing LT‐HSCs are modeled with a logistic population growth model (see Methods and Supporting Information). The total number of LT‐HSCs is believed to be fixed over the lifetime of an individual, and only a fraction are in the cell cycle at any time. The carrying capacity, N, in our model is interpreted as the homeostatic number of LT‐HSCs that are in the cell cycle (i.e., capable of cell division). These LT‐HSCs differentiate into short‐term hematopoietic stem cells (ST‐HSCs) with a fixed first order rate.

Each downstream compartment from ST‐HSC through colony‐forming unit erythroids (CFU‐Es) amplifies the influx of cells entering the next compartment (Figure 1a). These mechanisms are implemented using a system of ODEs described in detail in Methods and Supporting Information. Briefly, each compartment is characterized by two parameters: an amplification, , and a mean residence time, τ (Figure 1a and Table 1). The amplification parameter controls the degree of cell proliferation within a compartment. At steady‐state, the rate of cells leaving a compartment is amplification times the rate of cells entering that compartment (efflux = a × influx). The mean residence time (τ) is approximately how long a cell spends in that differentiation stage (see Supporting Information for a more precise mathematical definition). The ratio log2(a)/τ is approximately the mean rate of cell division within a compartment.

Hemoglobin production was modeled by an embedded set of molecular reactions within the RET and RBC compartments (Figure 1b). We considered four major hemoglobin types: HbA and HbS, the normal and sickled forms of the dominant hemoglobin in adults, fetal hemoglobin, HbF, and HbA2 (α2δ2). Each hemoglobin consists of four globin monomeric subunits. The model implements de novo monomer synthesis at a constant rate within reticulocytes. Monomers assemble into heterodimers and then hemoglobin tetramers, as shown in Figure 1b. No new monomers are synthesized within RBCs, but hemoglobin assembly reactions continue. Dimers and tetramers are assumed to have the same lifetime as RBCs. All reaction rate constants are listed in Table 2.

In the body, RBC production is regulated by a feedback loop via erythropoietin (epo). Epo is synthesized by the kidneys and liver in response to hypoxia, and targets CFU‐E to stimulate RBC production. , , To model this feedback, the model uses hemoglobin levels to compute venous blood oxygen, and maps this onto the CFU‐E amplification. This empirical mapping was parameterized to match the estimated rates of RBC production in healthy and SCD individuals (Figure 1c, and Methods).

Simulated healthy and sickle cell disease phenotypes

Model parameterization for healthy and SCD phenotypes was established using a combination of literature review, parameter tuning, and sensitivity analyses (see Tables 1 and 2, and Supporting Information for details). The key difference between the two phenotypes is the shorter RBC lifespan in SCD (12 days) compared to healthy (120 days), and an elevated RBC production (≈150 times amplification at the CFU‐E stage in SCD versus ≈31 times in healthy erythropoiesis at steady‐state) due to the epo‐driven feedback described above.

The system was initialized with a single LT‐HSC and evolved to steady‐state for each parameterization (Figure 4). At steady‐state, the healthy phenotype contains 24 trillion RBCs and ~13 g/dL of hemoglobin, whereas the SCD phenotype has half as many RBC (12 trillion) and ~7 g/dL of hemoglobin (Figure 4 top and middle panels). These values are within the reference clinical range for each condition. , The model computes venous blood oxygen for each phenotype using the predicted hemoglobin profile (Figure 4 bottom panels, see Methods for details). Steady‐state vO2 levels are within the reported reference range for each phenotype. The SCD steady‐state was used as the pretreatment baseline for treatment simulations described next.

FIGURE 4

Simulated erythropoiesis to steady state for healthy (left) and SCD (right) phenotypes using the endogenous model. Each simulation was initiated with a single LT‐HSC and no other cells. Simulation outputs shown are cell counts (top panel, 1 K = 103, 1 M = 106, 1B = 109 and 1 T = 1012), hemoglobin levels (middle panel), and venous blood oxygen level (bottom panel). The steady state value for each readout is shown in the plot legends. BFU‐E, burst‐forming unit‐erythroid; CFU‐E, colony‐forming unit erythroid; LT‐HSC, long‐term hematopoietic stem cell; MPP, multipotent progenitor cell; RBC, red blood cell; RET, reticulocyte; SCD, sickle cell disease; ST‐HSC, short‐term hematopoietic stem cell

Model of autologous stem cell transplantation gene therapy

We expanded the model to simulate a gene therapy for SCD using autologous stem cell transplantation. The therapy isolates CD34+ cells from a patient for ex vivo transduction with a lentiviral vector, and transplants treated cells back into the patient. Figure 2a schematically shows our implementation of the initial transplantation and subsequent engraftment of the modified progenitor cells.

CD34+ cells comprise multiple progenitor cell types. Each of the compartments from LT‐HSC through common myeloid progenitor (CMP) is seeded with a fraction of the total transplanted pool of CD34+ cells as shown in Figure 2a (also see Table 3). Each progenitor pool of transduced cells is assigned a separate branch. With this setup, treatment simulations can trace the time evolution of each branch and measure the contribution of that precursor pool to peripheral RBCs over time.

The combined pool of endogenous, non‐transduced, and transduced LT‐HSC follows the same logistic growth model as the pretreatment endogenous LT‐HSC pool. This means that the growth rate of the engrafted LT‐HSC is assumed to be identical to that of the endogenous pool, and the total available niche for all LT‐HSC is unchanged. The model also treats endogenous LT‐HSC and non‐transduced LT‐HSC in the drug dose identically.

Cells in the transduced branches contain the treatment gene. As a prototypical example, we considered an exogenous γ‐globin gene that encodes a novel γ‐globin (γNew) downstream of an endogenous β‐globin promoter. γNew is synthesized in RET derived from transduced precursors. This initiates an additional series of molecular reactions producing the novel fetal hemoglobin HbNew (Figure 2b). These reactions are analogous to model reactions for endogenous hemoglobin assembly shown in Figure 1b.

Before transplantation, a patient’s endogenous precursor cells are depleted with a bone marrow preconditioning regimen. The objective of this myeloablative preconditioning is to create a sufficient niche for the engraftment of transplanted cells. We approximated preconditioning by depleting the LT‐HSC, ST‐HSC, CMP, and multipotent progenitor cell (MPP) compartments to a fraction of their pretreatment steady‐state levels (Figure 2c, Table 3). Simultaneously, the system is seeded with the initial dose of ex vivo transduced autologous CD34+ cells.

Simulated treatment

Treatment model parameterization was informed and validated by published literature, as described in detail in the Supporting Information. Simulated post‐treatment time course for a patient with nominal SCD after applying preconditioning treatment and transplantation regimens described above is shown in Figure 2d,e. Treatment efficacy is measured by transduced RBC numbers, hemoglobin levels, and venous blood oxygen (Figure 2d). Simulations show a significant depletion of the endogenous RBC pool within the first few weeks after the transplantation due to the loss of endogenous precursors after preconditioning (Figure 2c). Over the next 1 to 3 months, as the progenitor pools are replenished, both endogenous and transduced RBC begin to appear in the blood. Over 6 to 12 months post‐transplantation, a new steady‐state is established with an overall higher number of RBC relative to the pretreatment baseline (top panel in Figure 2d, also see bottom panel of Figure 2e). The higher post‐treatment steady‐state population reflects the longer lifespan of transduced RBCs.

Endogenous hemoglobin levels drop with RBC numbers in the first few weeks post‐treatment (a patient may typically be given compensatory blood transfusions during this period). As the RBC numbers recover, so does the total hemoglobin (Figure 2d, middle panel). The post‐treatment dynamics of the therapeutic hemoglobin, HbNew (orange curve) show an early, transient phase in the first 3–6 months post‐transfusion, followed by a slower, sustained phase that approaches steady‐state in approximately a year. As explained below, this biphasic behavior is due to two waves of RBC derived from short‐lived and long‐lived progenitors.

The total hemoglobin level at the post‐treatment steady‐state is higher than the pretreatment baseline and consists of ~ 30% of the therapeutic hemoglobin (HbNew). As a functional effect of this, the model predicts an elevated blood oxygen level of 8.7 ml/dl at the post‐treatment steady‐state compared to 5.7 ml/dl at the pretreatment baseline. The increase in blood oxygenation is the result of both the increase in total hemoglobin, and the higher oxygen carrying capacity of the fetal isoforms HbF and HbNew.

Short and long‐term dynamics

To resolve the short and long‐term post‐treatment dynamics, we examined the engraftment time course of different progenitor pools (Figure 2e, top panel). The model assumes that vector‐transduced LT‐HSC have the same replication and differentiation rates as endogenous LT‐HSC, and these cells compete for the same bone marrow niche. In other words, the bone marrow microenvironment does not distinguish between endogenous and transplanted LT‐HSC. As a result, the fraction of transduced LT‐HSC is established at the time of the transplantation, and stays constant thereafter (Figure S2). The total number of LT‐HSC at steady‐state approaches the set carrying capacity of that compartment in the model.

In contrast to the LT‐HSC, short‐lived progenitor cells (ST‐HSC, MPP, and CMP) in the initial transfusion are completely depleted within the first 6 months. This is the expected model behavior, as these compartments have a limited regenerative capacity (i.e., a fixed number of cell divisions), and must be replenished by an upstream source of cells—ultimately, the LT‐HSC pool.

The transient engraftment of the short‐lived progenitors (ST‐HSC, MPP, and CMP) and stable engraftment of LT‐HSC in the initial dose produce two waves of transduced RBCs in the blood (middle and bottom panel of Figure 2e). The first wave, driven by short‐lived precursors, appears within 1 to 3 months post‐transplantation, and then tapers off over the next 9–12 months as these precursors terminally differentiate and die off. The second wave, driven by the LT‐HSC, is slower to arrive, but over 6 to 12 months post‐transplantation it establishes stable, long‐term erythropoiesis producing transduced RBCs. This biphasic behavior has been described in the literature. , In our model, it naturally arises from the self‐regenerating property of the LT‐HSC compartment in contrast to other progenitor pools that are eventually depleted without an upstream source.

Sensitivity to treatment parameters

To assess how changes in treatment parameters affect model predictions, we performed several sensitivity analyses (Figure 3 and Figure S3). The impact of treatment parameters most likely to affect key long‐term clinical outcomes is shown in Figure 3. Increasing the initial cell dose, the proportion of LT‐HSC in the dose and the mean vector copy number (VCN) all positively affect the long‐term outcome, as expected. Depleting greater than 90% of the endogenous LT‐HSC in the preconditioning stage is also predicted to increase treatment efficacy by opening up a bigger niche for engraftment. Thus, the intensity of the initial myeloablation is likely to be a critical determinant of therapeutic efficacy in the clinic (Figure 3, preconditioning panels). The effect of varying the proportion of other progenitors on short‐term outcomes is shown in Figure S3. The greatest impact is due to a higher proportion of ST‐HSCs.

We also used the sensitivity analysis to understand how RBC survival times affect model predictions. RBCs derived from transduced progenitors likely contain variable amounts of sickling HbS and anti‐sickling HbF (HbNew). Therefore, the lifespan of these RBCs will likely be shorter than the 120 days for normal RBCs expressing exclusively HbA. The sensitivity analysis predicts that increasing lifespan leads to greater levels of HbNew, fraction of non‐sickled RBC and venous oxygen in peripheral blood (Figure 3, right panel). Clinical studies have shown that a threshold of 10% HbF reduces major organ damage, 20% HbF results in amelioration of vaso‐occlusive crises and pulmonary complications in SCD, , and a 70% non‐HbS RBC alleviates the impact of sickle RBCs. Preclinical studies with anti‐sickling globins have shown complete correction of the SCD phenotype owing to a large selective advantage of HbF‐containing sickle erythrocytes over uncorrected ones, and a 20%–40% gene‐modified HSC chimerism results in a near pan‐cellular F‐cell production. Therefore, a lifespan of greater than or equal to 90 days of gene‐modified RBCs is expected to result in significant clinical improvement. However, future nonclinical data will be needed to adjust our model based on observed lifespan or survival rates of gene therapy treated RBCs.

DISCUSSION

In this study, we developed a multiscale mathematical model of erythropoiesis and gene therapy with autologous stem cell transplantation. We used the model to simulate the treatment of SCD with a novel exogenous γ‐globin that produces an anti‐sickling hemoglobin. SCD arises from a defect in a single gene, and is currently incurable in a large majority of patients. It is therefore an attractive target for gene therapy. The insertion of an alternative γ‐globin gene is shown to confer anti‐sickling behavior in preclinical studies. However, there is limited clinical data to provide a meaningful quantitative understanding of how treatment parameters influence post‐treatment outcomes. Our model was aimed at understanding these relationships and creating a rational framework to design and optimize gene therapies for SCD.

The model uses ODEs to implement three sets of interlinked mechanisms: (i) RBC production, (ii) hemoglobin synthesis, and (iii) regulation of RBC production by blood oxygen. Steady‐state behavior of healthy and SCD phenotypes in the model is consistent with clinical observations. The treatment model implements bone marrow preconditioning, stem cell engraftment, induction of erythropoiesis from transduced stem cells, and the expression of the anti‐sickling hemoglobin in peripheral blood. Treatment simulations resolve the early phase of reconstitution driven by short‐term progenitors in the first 6 months post‐engraftment, and the establishment of stable erythropoiesis driven by long term HSC over 12–18 months post‐engraftment, as reported in the literature.

Our approach provides a conceptual framework to understand and optimize gene therapy for SCD. For example, we used a sensitivity analysis to understand how initial cell dose and transduction efficiency affect long‐term treatment outcome. Successful treatment depends on efficient engraftment of treated LT‐HSC that can provide a stable, long term source of functional RBC, whereas a poor engraftment leads to reduced platelet and neutrophil recovery and negatively impacts patient survival. , The recommended minimal dose of stem cell transplantation for clinical studies is ≈2–5 × 106 viable CD34+ HSPCs/kg body weight, as this threshold is believed to produce sufficient engraftment. , Nominal treatment simulations with our model assumed a total dose of 108 cells (= 2 × 106 cells/kg × 50 kg). A sensitivity analysis of the model around this value predicts that increasing dose beyond this threshold will continue to increase the proportion of anti‐sickling hemoglobin.

However, substantially increasing the dose is limited by the challenge of harvesting sufficient CD34+ cells for ex vivo transduction. A sensitivity analysis with respect to VCN predicts that with the nominal dose fixed, increasing transduction efficiency (i.e., mean VCN in our model), is also likely to produce a similarly beneficial long‐term outcome as a higher dose. The analysis also reveals the critical role of the myeloablative preconditioning to create a sufficient niche for successful engraftment. These sensitivity analyses thus permit a meaningful cost‐to‐benefit analysis to optimize the combination of dose, transduction efficiency, and preconditioning to achieve a desired outcome. Model readouts, such as the proportion of transduced RBCs, can be related to clinical biomarkers such as the proportion of F‐cells. Model parameterization and predictions can therefore be further refined with emerging data from early clinical trials, and help guide treatment optimization.

To limit model complexity, we used several simplifying assumptions. First, we limit the vector copy number per cell to one, and treat mean VCN as the fraction of cells containing the vector. This assumption ignores the heterogeneity of vector distribution, and, in particular, it ignores the possibility of a mean VCN greater than 1 due to cells taking up multiple copies of the vector during ex vivo transduction. Second, we model the effect of preconditioning prior to treatment transfusion by depleting various progenitor pools by a fixed fractional amount. This simplification ignores the pharmacokinetics of the preconditioning drug, and the variable depletion of progenitor cell types. All infused cells are assumed to home rapidly to the bone marrow. As the focus of this study is SCD, which affects RBCs, we ignored other branches of hematopoiesis and the production of other cell types in blood. Thus, the current model does not account for any changes in the relative production rates of different cell types upon an acute perturbation from steady‐state hematopoiesis. Another simplification was to ignore the association of heme groups with globins, and the role of iron in hemoglobin formation. Finally, the feedback regulation of RBC production via blood oxygen was implemented empirically without a specific mechanism. These simplifications allowed us to focus on the most significant pieces of a complex biological process, and represent them within a deterministic, mathematical structure.

The QSP model presented here is intended to be a predictive tool for optimizing stem cell transplantation‐based gene therapies for SCD and other genetic disorders. We demonstrated how to represent the complex biology of erythropoiesis and hemoglobin production within a multiscale ODE model. We used this system to understand how gene therapy treatment parameters are quantitatively linked to long‐term outcomes. Several simplifications used in the current approach are noted above. In future studies, we intend to refine the model further to address these.

CONFLICT OF INTEREST

L.W., D.H., A.M., J.S., F.H., J.F.A., J.M.B., and R.D. are employees of Applied BioMath LLC. B.Z., K.P., J.A., T.Y., E.C., and J.R., are employees of CSL Behring, King of Prussia, PA, USA. All other authors declared no other competing interests for this work.

AUTHOR CONTRIBUTIONS

B.Z., L.W., and R.D. wrote the manuscript. B.Z., L.W., and R.D. performed the research. L.W. and R.D. analyzed the data. B.Z., L.W., K.P., D.H., A.M., J.A., J.S., F.H., T.Y., E.C., J.F.A., J.M.B., J.R., and R.D. designed the research.

Supplementary Material

Click here for additional data file.