Modeling changes in biomarkers in Gaucher disease patients receiving enzyme replacement therapy using a pathophysiological model.

BACKGROUND
Gaucher disease (GD) is a rare recessively inherited disorder caused by deficiency of a lysosomal enzyme, glucocerebrosidase. Accumulation of glucosylceramide or glucosylsphingosine in macrophages leads to increased production of ferritin and chitotriosidase and to decreases in hemoglobin concentration and platelet count, which are used as blood biomarkers. GD is treated by enzyme replacement therapy (ERT) or, sometimes by substrate reduction therapy. However, no physiological model for analysis of biomarkers change during ERT has been proposed. We aimed to develop a pathophysiological model to analyze biomarker's response to ERT and several covariates impact.


METHODS
Changes in blood ferritin, chitotriosidase, hemoglobin and platelets were analyzed in French GD Registry patients receiving imiglucerase/alglucerase as ERT. We used simplified exponential pathophysiological model, with initial concentration, biomarkers amplitude of variation and rate constant of normalization during ERT. Changes in four biomarkers were analyzed separately and then all four together from initiation to discontinuation of ERT, or until the end of follow-up. Several covariates were tested, including age at ERT initiation, splenectomy, sex, genotype (N370S/N370S), and ERT dose.


RESULTS
An exponential model gave a good data fit. The four biomarkers analysis showed that the rate of nomalization was the same for all biomarkers, with a half-life of 0.5 years. Predicted values of biomarkers at ERT's steady state were 40% and 10% of initial concentrations, for ferritin and chitotriosidase, respectively, and 120% and 200% for hemoglobin and platelets, respectively. We found that 3 covariates had an effect on initial concentration or on amplitude of variation in ferritin, hemoglobin and platelets: women and patients under 15 years of age had lower ferritin and hemoglobin concentrations, and patients under 15 years of age had higher platelet count. Splenectomized patients had higher ferritin concentrations and platelet count and lower amplitude of variation of hemoglobin.


CONCLUSION
We report the first dynamic model of biomarker changes in GD. It enabled us to estimate that 95% of biomarker response to ERT was achieved in 2 years, but with high inter-patient variability. We also found that with the current treatment, normalization of chitotriosidase and ferritin will occur in about 65% of patients.


Background
Gaucher disease (GD) [1] is a recessively inherited lysosomal storage disorder caused by deficiency of a lysosomal enzyme, glucocerebrosidase (EC 3.2.1.45), which leads to insufficient clearance of the enzyme's substrate, cellular glucosylceramide. Pathologic accumulations of glucosylceramide (or other substrates, such as glucosylsphingosine) in the lysosomes of tissue macrophages (Gaucher cells) results in splenomegaly, hepatomegaly and multiple forms of skeletal disease [2]. Three clinical phenotypes have been described: type 1, the prevalent form usually defined by the absence of central nervous system impairment; and types 2 and 3, both rare and severe, have central neurological involvement [3]. GD diagnosis is confirmed by the detection of low glucocerebrosidase activity, usually less than 30% of the normal value in peripheral leukocytes. Genotyping can sometimes provide prognostic information [4]. More than 250 mutations of the GBA1 gene encoding lysosomal glucocerebrosidase have been reported as being associated with GD, but the predominant mutation in type 1 GD is called N370S (or c.1226A > G) [5]. The N370S mutation is usually protective against neuronopathic disease. GD can be treated by enzyme replacement therapy (ERT). The first enzyme preparation used to treat GD consisted of placenta-derived glucocerebrosidase (alglucerase available in 1991, Genzyme Corporation) with modified mannoseterminated glycans, allowing more selective uptake by tissue macrophages, the prominent storage cells in GD [6][7][8]. This preparation was replaced in 1996 by recombinant enzyme (imiglucerase), which was therapeutically equivalent in terms both of safety and efficacy [9]. Two new biosimilar agents are now available: velaglucerase-alfa (Shire) [10] and taliglucerase-alfa (Pfizer) [11]. ERT reduces macrophagic substrate accumulation, but no routine substrate assay is currently available. A substrate reduction therapy (miglustat, Actelion) that can be prescribed in a few indications has been available since 2002.
The levels of several biomarkers (e.g., in our article, ferritin and chitotriosidase, but also tartrate-resistant acid phosphatase or angiotensin-converting enzyme, not analyzed in this study) change during the clinical course of GD [12][13][14] due to macrophagic activation: their concentrations rise with disease progression and generally decrease during ERT [12,15,16]. These variations can predict bone complications [17]. High blood ferritin during the course of type 1 GD may reflect macrophage activation triggered by substrate accumulation, demonstrated by increases in CCL18 and macrophage inflammatory protein-1α or 1β [18,19]. ERT is associated with a dramatic decrease of blood ferritin [16,17], which is more pronounced in patients with an intact spleen [20]. Chitotriosidase is massively produced by storage cells and there is a linear relationship between chitotriosidase and glucosylceramide levels, as shown in spleen sections from patients with GD [21]. Chitotriosidase values drop sharply during ERT, when substrate accumulation decreases, coinciding with clinical improvements [12].
Hematological abnormalities (anemia and thrombocytopenia) are common in GD, because Gaucher cell infiltration leads to hypersplenism (increased destruction or sequestration of red blood cells or platelets) and bonemarrow insufficiency (decreased production) [22]. Splenectomy, performed essentially before 1991, increases baseline platelet count and decreases the slope of platelet clearance during ERT [17]. Anemia and thrombocytopenia can be used as biomarkers to manage GD patients.
Grabowski et al. [23] developed an Emax model to describe changes in hemoglobin and platelets and in splenic and hepatic volume during ERT of patients in the International Collaborative Gaucher Group Registry. Biomarkers of French GD patients from a single center was modeled before and during ERT by Stirnemann et al. [17], but no physiological model was proposed to analyze changes in biomarkers levels during ERT. Nonlinear mixed effects models [24] are widely used to analyze biological processes described by repeated longitudinal data. They allow estimation of the mean value of the parameters and their interindividual variability. These models allow a sparse sampling design with few data points per individual in a large set of individuals.
The aim of this study was to develop a pathophysiological model explaining the response of biomarkers (ferritin, chitotriosidase, hemoglobin, platelets) to ERT and to analyze the influence of several covariates.

Patients and data
We analyzed changes in four biomarkers, ferritin, chitotriosidase, hemoglobin and platelets, in patients receiving ERT with alglucerase from 1991 to 1996 and thereafter imiglucerase, from the French GD Registry (FGDR) [25]. The French Data Protection Commission (CNIL) approval of the FGDR required oral or written informed consent from patients or their parents. The local Institutional Review Board of Northern Paris Hospitals, Paris-Diderot University, AP-HP (Ethics Committee) reviewed and approved the initial research project (220-08). These patients were included in the study if the last follow-up visit was between 2009 and 2010. This investigation was undertaken to describe changes in biomarkers from initiation of ERT until 31 December 2010 corresponding to the limit of collected retrospective data in the FGDR. These data correspond also to those used in the first description of patients of the FGDR [25]. All measurements of biomarkers from 2 years before the initiation of ERT to discontinuation of ERT (interruption of more than 6 months) or the end of followup (December 2010) were used. When the chitotriosidase activity was undetectable (patients with probable homozygous chitotriosidase-gene deficiency), this biomarker was not retested [26] and the data were not included in the analysis. When ERT was interrupted for less than 6 months, patients were considered to be still under treatment. Several covariates were tested including age at initiation of ERT, splenectomy, genotype (N370S/N370S or others), sex, dose at initiation of treatment (divided into 3 classes: lower than 90 IU/kg/month, between 90 and 120 IU/kg/month or higher than 120 IU/kg/month) and average dose in the third year of treatment (divided into the same 3 classes) for patients followed up for at least 3 years.

Model for biomarker changes during ERT
We defined a pathophysiological model [27] of GD and treatment by ERT ( Figure 1). GD is caused by a glucocerebrosidase deficiency that leads to an accumulation of glucosylceramide, which is no longer degraded. As a consequence, there is an increase in ferritin and chitotriosidase production. With ERT, patients are supplied with glucocerebrosidase which leads to the degradation of the glucosylceramide and then improvements of biomarkers levels. The differential equation of this model for chitotriosidase is presented in Supplementary material (Additional file 1). As only biomarker measurements were available, we made a further simplification. Since the rate constant of biomarker elimination is very rapid compared with the rate constant of normalization of glucosylceramide under treatment (k), it was neglected so that we have an exponential increase for chitotriosidase. For chitotriosidase activity, the model is: where C 0 is the initial concentration, r is the amplitude of variation and k is the rate constant of normalization during ERT. For serum ferritin, a similar reasoning leads to the same equation. From the model, we can derive the corrected values of biomarkers at steady state of ERT as C 0 r ( Figure 1) and the normalization half-life as log 2 ð Þ k . Even though two different mechanisms explain hematological abnormalities (destruction by hypersplenism and decrease hematopoiesis), we proposed to simplify the model with only one predominant mechanism (destruction). By similar reasoning, for hemoglobin and platelets, the model is: which increases from C 0 to C 0 r ( Figure 1).

Statistical analyses
We used a nonlinear mixed effects model to analyze measurements of the four biomarkers. First, we performed a separate analysis of each biomarker. We assumed an biomarkers that increase during ERT (hemoglobin, platelets). C 0 is initial concentration, k rate of decrease (or increase) and C 0 r estimated corrected value at steady state of ERT. exponential random effect on all parameters and that these random effects have a normal distribution with a mean of 0 and a standard deviation of ω. The residual error model was supposed to be proportional with a standard deviation, σ. We also tested the correlation between random effects. Model selection was performed using the Bayesian information criterion (BIC) which adds to the log-likelihood penalty increases with the number of model parameters and the sample size used (BIC = −2 × logL + k log (N)), where L is the maximum likelihood, k is the number of model parameters and N is the number of patients. BIC is a parsimonious criterion where the lowest value corresponds to the best model. We compared our exponential model with the Emax model of Grabowski et al. [23] using the BIC.
We tested the impact of the covariates mentioned above to explain part of the variability [28]. As only discrete covariates were considered, the effect of the covariate is to change the parameter by exp(β) where β is the estimated regression coefficient. The final model was built using a two-step approach. In the first step, individual empirical Bayesian estimates of parameters were generated from the basic model without covariates. We searched for univariate associations with covariates using Wilcoxon or Spearman tests. In the second step, all covariates with a p-value of < 0.2 were entered into a multivariate model. We then performed a multivariate analysis using a forward selection. Pvalues of the Wald test were assessed at the 0.05 level. Splenectomy has an impact on changes in platelet count [17]. So we considered (and tested) that splenectomy has an impact on the initial platelet count and on its amplitude of variation.
We then analyzed the 4 biomarkers jointly and evaluated whether one or more parameters were correlated or if several parameters had similar values, to simplify the model. We evaluated the final model with various goodness-of-fit plots [29]: individual fits, graphs of individual weighted residuals versus time and visual predictive check. The significant covariates obtained previously in models of each biomarker were added. Then, we performed a backward selection and computed the P-values of the Wald test. Following estimation of population parameters, we could estimate individual parameters. Individual empirical Bayesian estimates were obtained as maximum a posteriori. From the individual parameters, we estimated individual half-life, time to 95% of response (as 4.3 half-lives), values of biomarkers at steady state of treatment and whether this value is normal or not.
Estimations were performed using the Stochastic Approximation Expectation Maximization (SAEM) algorithm in MONOLIX 4.2.0 (Lixoft, Orsay, France, available at http://www.lixoft.com) [30]. BIC and log-likelihood were estimated by importance sampling.  hemoglobin and platelets, respectively. Median length of follow-up during ERT was 5 years for ferritin and 6 years for the other biomarkers. Over 65% of the patients were followed up for more than 3 years. Figure 2 shows the individual changes in each biomarker during ERT in a spaghetti plot. The exponential model fitted the data. Comparison of our exponential model with the Emax model gave a gain BIC of 3 for ferritin, 12 for chitotriosidase, and 22 for platelets for our model and a loss of 7 for hemoglobin. We tested the impact of covariates in the model of separate biomarkers and found that 3 covariates (initiation of ERT < 15 years of age, splenectomy and sex) had an effect on the initial concentration and/or on the amplitude of variation in ferritin, hemoglobin and platelets. Genotype and dose had no significant impact on the parameters and none of the covariates had a significant impact on chitotriosidase. The best model of analysis of the four biomarkers together was the one with the same parameter k for all biomarkers, so we considered a rate constant of normalization during ERT common to the four biomarkers. The model with the same parameter k had a smaller BIC of 22 than the model with 4 different values for parameters k. Values of parameters without covariates (only splenectomy for platelets) are given in Table 2 and Table 3: Predicted values of biomarkers at steady state of ERT were 38% and 10% of initial concentrations, for ferritin and chitotriosidase, respectively, and 117% and 200% for hemoglobin and platelets, respectively. Estimated corrected values at steady state of ERT were 179.4 ng/L, 792 nmol/h.mL and 13.9 g/dL for ferritin, chitotriosidase and hemoglobin, respectively. The estimated corrected platelet counts at steady state of ERT were 156,030/mm 3 and 195,037/mm 3 for nonsplenectomized and splenectomized patients, respectively. Individual predicted values show that, under current ERT, normalization will occur for 58% of the patients for ferritin, 66% for chitotriosidase, 88% for hemoglobin and 64% for platelets. Figure 3 shows σ_P 20 3 C 0 is initial concentration, r amplitude of variation and k the rate of normalization during ERT. F corresponds to ferritin, C to chitotriosidase, H to hemoglobin and P to platelets. The P-value is for the Wald test. Exp β represents the effect of the covariates. Splen corresponds to the covariate splenectomy. ω: standard deviation of inter-individual variability, σ: standard deviation of proportional residual error. RSE: relative standard error. For example, baseline ferritin concentration is C 0F =598 μg/L, and the estimated corrected value at steady state of ERT is C 0F × r F = 179.4 μg/L. Baseline platelet count in non-splenectomized patients is C 0P = 74,300/mm 3 and C 0P exp(β _ splen) = 185,750 /mm 3 in splenectomized patients. The estimated corrected platelet count at steady state of ERT in splenectomized patients is (C 0P × exp(β _ splen)) × (r 0 × exp(β _ splen)) = 195,037 /mm 3 . individual fits for one splenectomized patient and one non-splenectomized. Relative standard errors were lower than 20% for both fixed effects and variabilities. Diagnostic plots (Additional file 2) show that the model describes the data adequately. Values of parameters in the joint model with covariates are given in Table 4. The covariate effects were the same as in the models considering the biomarkers separately. Figure 4 shows the response of the biomarkers to ERT, according to significant covariates. Non-splenectomized (p-value = 10 −6 ) and age <15 years (p-value < 10 −10 ) patients had a lower initial concentration of ferritin. Women (p-value = 0.002) and non-splenectomized (p-value = 10 −4 ) patients had a lower amplitude of variation in ferritin. Women (p-value = 0.01) and age <15 years (p-value = 10 −8 ) patients had a lower initial concentrations of hemoglobin. Splenectomized (p-value = 10 −5 ) patients had a lower amplitude of variation in hemoglobin. Splenectomized (p-value < 10 −10 ) and age <15 years (p-value = 10 −6 ) patients had a higher initial platelet count. Splenectomized (p-value = 10 −8 ) patients had a lower amplitude of variation in platelets.

Changes in biomarkers
In the final model, we estimated that during ERT the half-life of normalization was 0.5 years (95% CI = [0.4-0.6]), with a variability of 97%. Half-life was 0.5 (0.1-3.1) years and 95% of response was obtained after 2 (0.4-13.3) years under ERT. Variabilities were reasonable for ferritin, hemoglobin and platelet count. Greater variability was found for chitotriosidase, doubtless because of variability in measurement error (a sample dilution is usually needed when values are high).

Discussion
Our pathophysiological model predicts changes in biomarkers on ERT and estimates the rate constant of normalization. For the final model, we estimated a normalization half-life of 0.5 years during ERT for all four biomarkers. Only 2 studies have modeled the changes in biomarkers and they used different models: Emax and linear mixed models. In the study of Grabowski et al. [23], the ERT dose effect had a significant impact on biomarkers with a large sample size. Emax is an empirical model with a hyperbolic function; it is not the result of a physiological model. These models are used in pharmacometrics [27] to study the link between dose and concentration. However, the drug often acts on a biological quantity, modifying its production or elimination. We modeled this quantity by a model with a production rate and a rate constant of elimination. Using physiological knowledge, we obtained differential equations to explain changes in biomarkers over time. For instance, 90% response is obtained after 9 T 50 for an Emax model, whereas it is reached after only 3.3 half-lives (similar to T 50 ) for an exponential model.
Our results show improvement in all biomarkers under ERT (decrease in ferritin and chitotriosidase and increase in hemoglobin and platelets). Stein et al. [20] also highlighted an increase in ferritin and a renormalization under ERT. Hollak et al. [31] reported a decrease of chitotriosidase of 32% in 1 year, and we found a 95% response in 2 years and a 36% response in 1 year. De Fost et al. [32] reported that 53% of patients had anemia at baseline and 58% had thrombocytopenia, with renormalization under ERT, and showed a similar pattern of response after 1 year under ERT.
Patients <15 years of age have lower initial concentrations of ferritin and hemoglobin but higher platelet counts; at initiation of ERT, women have a lower concentration of hemoglobin and splenectomized patients have a higher platelet count and ferritin. For GD children from the International Collaborative Gaucher Group Registry, Kaplan P et al. [33] noted that 50% had platelet counts less than 120 × 10 3 /mm 3 and 40% had anemia at the time of diagnosis. Hemoglobin and ferritin tend to be lower in GD women, as for other women, probably because of menstruation, with no link to GD. In our model, no covariates had a significant impact on the chitotriosidase changes. In contrast, Stirnemann et al. [17] found that splenectomy and GD genotype affected chitotriosidase activity, but their study was limited by a small number of samples, and they used a different model.
Our model allows estimation of the rate constant of biomarker improvement using a pathophysiological model. Bone events are the most debilitating and disabling complication of GD. With substrate overload, Gaucher cells activate and induce proinflammatory cytokine synthesis which can modify the activity of the osteoblast-osteoclast system and promote lytic phenomena and intraosseous vascular complications [34,35]. Further analysis of the interaction between biomarkers and bone events is needed.
Recently, plasma glucosylsphingosine has been proposed as a biomarker for GD [36,37] and could be used as a reflection of intracellular glucosylceramide. Our model may also predict changes in intracellular glucosylceramide and/ or glucosylsphingosine.
Our model could be used to study the effect of biomarker changes on complications such as the occurrence of bone events. Using a few measurements of biomarkers post-treatment, we can estimate the rate constant of normalization and individual amplitude of variation in biomarkers in order to predict further response to treatment. Our model could help to define an individual risk for complications and to refine the best ERT regimens.
Parameters of this pathological model could be used to predict value at steady state with the dosages in the first months of ERT. These predicted values could be used as an individual objective of treatment to modify dosage for each patient. Clinical studies are needed to confirm this hypothesis.
We developed a model which could in the future be used to manage patients, but further studies are needed to confirm clinical applications. This model with 2 parameters (in addition to the baseline value) predicts the level of improvement of biomarkers and the time to 95% response, using only a few measurements per patient.

Conclusion
This is to our knowledge the first study of changes in biomarkers in Gaucher disease using a pathophysiological model. With the model we estimate that 95% of biomarker response to ERT is achieved in 2 years, but with high inter-patient variability. We also found that with the current treatment, normalization of chitotriosidase and ferritin will occur in about 65% of patients.

Additional files
Additional file 1: Equation describing change in chitotriosidase level.
Competing interests FM had consulting fees and INSERM UMR1137 received a grant from Sanofi. JS had travel fees from Sanofi-Genzyme. NB had travel fees, consulting fees, fees for speaking and received grants (Genzyme, Shire, Actelion, Pfizer) donated to the Department of Clinical Research of the Assistance-Publique Hôpitaux de Paris. The other authors declared no financial disclosures.
Authors' contributions MV, JS and FM designed the research and analyzed and interpreted the data. MV, JS and FM wrote the first draft of the paper, which was then corrected and approved by all authors.
(See figure on previous page.) Figure 4 Biomarker responses for a median patient for the final model according to significant covariates for ferritin, chitotriosidase, hemoglobin and platelets. A green line is for non-splenectomized patients and a red line for splenectomized patients. A dotted line is for man and continuous for woman. For platelets, no distinction between man and woman and change is a continuous line. Initial concentration and concentrations at steady state of ERT are indicated, in green for non-splenectomized, red for splenectomized patients, and in orange when no distinction can be made depending on splenectomy. Lines are thin for man and bold for woman. Left: patients <15 years of age. Right: patients > 15 years of age.