Modelling primaquineinduced haemolysis in G6PD deficiency
Abstract
Primaquine is the only drug available to prevent relapse in vivax malaria. The main adverse effect of primaquine is erythrocyte age and dosedependent acute haemolytic anaemia in individuals with glucose6phosphate dehydrogenase deficiency (G6PDd). As testing for G6PDd is often unavailable, this limits the use of primaquine for radical cure. A compartmental model of the dynamics of red blood cell production and destruction was designed to characterise primaquineinduced haemolysis using a holistic Bayesian analysis of all published data and was used to predict a safer alternative to the currently recommended once weekly 0.75 mg/kg regimen for G6PDd. The model suggests that a stepwise increase in daily administered primaquine dose would be relatively safe in G6PDd. If this is confirmed, then were this regimen to be recommended for radical cure patients would not require testing for G6PDd in areas where G6PDd Viangchan or milder variants are prevalent.
https://doi.org/10.7554/eLife.23061.001eLife digest
Malaria is the most important parasitic disease that affects humans. Over half of the malaria cases in Asia and South America are caused by a species of malaria parasite called Plasmodium vivax (known as vivax malaria). This form of malaria results in repeated illness because dormant parasites in the liver wake at intervals to infect the blood. The only available drug that can stop these relapses is a drug called primaquine, which was developed seventy years ago.
Unfortunately, primaquine causes dangerous side effects in certain individuals who are deficient in an enzyme called G6PD, which helps defend red blood cells against stresses. Primaquine damages these cells so that they burst, leading to anaemia. This is a major problem because G6PD deficiency is common in regions where malaria is present: in some areas up to 30% of the population may be G6PD deficient. Since G6PD testing is not widely available, doctors often avoid prescribing primaquine to treat malaria, which results in more cases of disease relapse. Failing to prevent vivax relapses causes extensive illness and hinders efforts to eliminate malaria.
Is there a way to give this drug to patients that would be safer for people with G6PD deficiency? Primaquine destroys older rather than younger red blood cells. Watson et al. used mathematical modelling to see whether it is possible to develop a primaquine treatment strategy that would allow a gradual destruction of older red blood cells in individuals with G6PD deficiency, which would be safer. The mathematical model incorporates data from previous studies in malaria patients and healthy volunteers with G6PD deficiency and combines this with knowledge of how red blood cells are produced and destroyed. Watson et al. predicted that giving primaquine over 20 days in a steadily increasing dose was safer than current recommendations.
Mathematical models are simplifications of real world processes. The only way to test these findings properly will be to run a clinical trial that gives healthy volunteers who are G6PD deficient a course of primaquine treatment with a steadily increasing dose.
https://doi.org/10.7554/eLife.23061.002Introduction
Radical cure of vivax malaria in G6PD deficient patients
Plasmodium vivax accounts for over half the world’s malaria burden outside subSaharan Africa (Gething et al., 2012). The control and elimination of vivax malaria require both cure of the blood stage infection (the stage that causes acute illness) and the prevention of later relapses which derive from dormant hypnozoites in the liver (radical cure). Hypnozoites are formed from sporozoites, which do not develop immediately following mosquito inoculation but instead remain dormant in hepatocytes for weeks or months before developing and causing recurrent blood stage infections called relapses. In general, P. vivax infections in tropical regions are associated with frequent relapses (with intervals as short as three weeks) whilst relapses in P. vivax infections from Central America, Northern India and temperate regions are associated with longer intervals from acute infection to first relapse (White, 2011).
Primaquine, an 8aminoquinoline, is currently the only widely available antimalarial drug for the radical cure of P. vivax infections. Primaquine causes predictable oxidant haemolysis in G6PD deficiency (G6PDd) one of the most common genetic abnormalities of man (Cappellini and Fiorelli, 2008). Throughout Asia, the Mediterranean littoral and Africa, allele frequencies for this enzyme deficiency vary between 3% and 35% in the populations at risk from vivax malaria (Howes et al., 2013). As G6PDd has sexlinked inheritance, males are either deficient (hemizygotes) or normal, whereas women can be deficient (homozygotes), normal or partially deficient (heterozygotes) in proportions determined by the HardyWeinberg equilibrium. Because of Lyonisation, there is substantial variability in the proportion of red cells which are deficient in individual heterozygote females (Beutler et al., 1962).
The degree of haemolysis following primaquine depends on the dose administered and the severity of the enzyme deficiency (and in heterozygote females the proportion of erythrocytes which are deficient). The more severe G6PDd variants found in South East (SE) Asia (for example, Viangchan, Mahidol, Coimbra, Union) and the Middle East/West Asia (for example, Mediterranean) are generally associated with more severe haemolysis compared to the common African A variant. For G6PD normal patients, the primaquine regimen for radical cure that is recommended in SE Asia and Oceania (where relapse rates are high) is 0.5 mg base/kg/day for 14 days. Elsewhere it is 0.25 mg/kg/day for 14 days. For patients with G6PDd, a weekly dose is recommended; 0.75 mg/kg/week given for a total of 8 doses. Unfortunately G6PDd testing is not widely available despite the recent introduction of pointofcare rapid diagnostic tests (RDTs) for G6PDd. These RDTs are currently too expensive to deploy on a wide scale and can be difficult to interpret, and thus are not generally available (Brito et al., 201608; Satyagraha et al., 2016; Oo et al., 2016). Thus, primaquine is commonly not given to patients to avoid the risk of haemolysis so the burden of vivax malaria remains high, causing considerable morbidity and economic loss (Price et al., 2007).
Mechanisms of red blood cell production
The mechanisms regulating red blood cell production and turnover have been well characterised. Red blood cells (RBCs) transport oxygen which is reversibly bound to the main red cell protein, haemoglobin. RBC production in the bone marrow is regulated to maintain oxygen carrying capacity. When the haemoglobin concentration in the blood falls, this reduces oxygen carriage and RBC production is upregulated, a process mediated largely by the renal hormone, erythropoietin. At times of increased bone marrow production, reticulocytes appear in increased numbers in the circulation (the upper limit of normal is $\approx 1.5\%$). Normal RBCs in healthy people have a very stable life expectancy of around 120 days. This is well modelled by a Gumbel distribution with low variance. In nucleated cells G6PD can be newly synthesised, but the red cells lose their nucleus before leaving the bone marrow so very young red cells (reticulocytes) have the highest G6PD activity, and this declines as the RBCs age. In most G6PDd variants, the mutant enzyme degrades more rapidly compared to the normal enzyme. Older erythrocytes may have up to five times less G6PD activity than reticulocytes. G6PDd results in lowered NADPH and a reduced ability to regenerate reduced glutathione. Reduced glutathione protects normal RBCs against oxidant stresses such as the haemolytic effects of primaquine metabolites and certain foods, classically fava beans. G6PD is also important for the function of catalase, another oxidant defence mechanism. As these nonreusable oxidant defence reserves are ‘used up’, the aging erythrocyte becomes increasingly vulnerable to oxidant haemolysis (Beutler et al., 1954a; Dern et al., 1954; Beutler, 2008; Recht et al., 2014).
Evidence from previous studies of oxidant haemolysis in G6PD deficiency
As young red cells have more functional enzyme than older cells, the degree of oxidant haemolysis depends on the genetic variant of G6PDd and the age distribution of the red cell population. Once the older cells have haemolysed, the remaining younger erythrocytes are essentially resistant to further damage by the same dosing regimen (that is, drug exposure) (Beutler et al., 1954a). However, higher primaquine doses do cause further haemolysis. This explains the fall then rise in haemoglobin with continued daily primaquine administration in mild and moderate severity variants of G6PDd. This temporary primaquine insensitivity in G6PDd individuals with the continued primaquine administration was characterised by Beutler and colleagues in a series of studies conducted over sixty years ago (Beutler et al., 1954a,1954b, 1955; Dern et al., 1954; Beutler, 1959) and later exploited by Alving et al. to develop the once weekly regimen in G6PDd (Alving et al., 1960.)
By experimenting with highdose weekly regimens and lowdose daily regimens, Beutler and colleagues showed haemoglobin would first fall as a result of oxidant haemolysis and then rise despite continued exposure to the same doses of primaquine which had caused the initial haemolysis. This resulted from reactive erythropoiesis (reticulocytosis) that introduced a younger red cell population to the circulation which was essentially ‘resistant’ to the haemolytic effects of that primaquine dose. Intermittent primaquine administration resulted in progressively smaller cycles of haemolysis followed by reticulocytosis as the red cell population became younger. These results led to a recommendation for a highdose, once weekly primaquine regimen for radical cure in vivax malaria patients with G6PDd (8 once weekly adult doses of 45 mg) (Alving et al., 1960). This regimen was devised based on studies in subjects with the African ${A}^{}$ variant of G6PDd, which is one of the mildest deficiencies. Safety was not formally assessed in more severe deficiencies. A recent trial of this regimen in vivax malaria patients with the more severe Viangchan G6PDd variant from Cambodia showed a greater fall in haemoglobin and a delayed recovery from anaemia in G6PDd compared to G6PD normal patients with one patient requiring a blood transfusion (Kheng et al., 2015). These data suggest that weekly primaquine may not be the optimal regimen for the more severe G6PDd variants prevalent outside Africa.
Reconsideration of the detailed haematological studies that laid the foundation for the weekly regimen suggests that an ascendingdose regimen of primaquine, with a schedule that matches the dynamics of red blood cell production, could induce a safe ‘slow burn’ haemolysis, even in individuals with severe G6PDd variants, and would still deliver a total therapeutic dose for radical cure.
Accordingly, our study had two objectives; first, to construct a compartmental model for red blood cell dynamics which could be used to analyse all available data from past studies of haemolysis in G6PDd individuals, and second to predict an optimal ascending dose regimen which would be safe and efficacious yet practical and could, therefore, be recommended without G6PD testing.
Results
Model fit
Figure 1 shows hypothetical data simulated from the compartmental model with a primaquine regimen of 45 mg weekly for eight weeks fitted to data from adult G6PD deficient Cambodian patients. Parameters were randomly drawn from the Bayesian posterior distribution.
The signaltonoise ratio in the reticulocyte data is low and this is apparent from the median reticulocyte count which varies considerably during the 56 days. In comparison, simulations from the mechanistic model show that a substantial rise in the reticulocyte count should occur approximately one week after the first dose, with a peak after the third dose, and then return to normal slowly over the subsequent six weeks. The serial haemoglobin data on the other hand show a clear trend with a large fall after the first dose, a smaller fall after the second and then a gentle recovery with no major effect from subsequent primaquine doses. This trend is reproduced by the model and the posterior distribution also characterises satisfactorily the variance observed in steady state haemoglobin concentrations.
Predicted dose response
Combining the data from Figures 2–4, it is possible to estimate a primaquine dosehaemoglobin response curve for G6PDd individuals whose severity is similar to the ‘moderate severity’ variants G6PDd Mahidol/Viangchan. The data at different dosing levels are sparse and the studies have been done in very different contexts; however, the strong mechanistic assumptions used to construct the compartmental model regularize the problem enough to compare the studies in a principled way. The data from G6PDd Mediterranean are excluded from this doseresponse curve estimation because the haemolysis observed with this variant is considerably greater than for G6PDd Mahidol/Viangchan. However, the observed falls in haemoglobin after 5 daily doses of 30 mg in G6PDd Med Sardinians are shown by the red triangles in Figure 5, right plot, for comparison.

Figure 4—source data 1
 https://doi.org/10.7554/eLife.23061.007
The posterior MCMC samples inferred from the Kheng data can be used to approximate model uncertainty around the median doseresponse curve. The right plot of Figure 5 shows the posterior predictive doseresponse curve with 90% credible intervals, where the ‘response’ is defined as the drop in haemoglobin after five days at a given dosing level. Overlaid are estimates of the falls in haemoglobin induced by 5 daily doses from studies in Figures 2,3, and an extrapolated estimate from the posterior distribution of the model fitted to data from weekly dosing in Viangchan variant.
It is of interest to compare the fitted doseresponse relationship in Figure 5 (right: thick black line)—corresponding to the more severe variants of G6PDd—with the green crosses corresponding to observed and fitted haemolysis in G6PDd African ${A}^{}$ (mild variant). As would be expected, for the mild variant the doseresponse relationship has the same shape but is shifted to the right.
Safe optimal regimen
The currently recommended dose for the radical cure of vivax malaria in an adult in SE Asia and Oceania delivers 420 mg (that is, 30 mg/d x 14 d) of primaquine and is very effective (John et al., 2012). The maximum primaquine dose administered in the weekly regimen is 360 mg (8 x 45 mg) but the efficacy of this regimen has only been reported in Afghan refugees in Pakistan, a country with a relatively low relapse rate (Leslie et al., 2008).
The primary objective of our research is to design a novel primaquine regimen that could be given safely to individuals with G6PDd or of unknown status without G6PD testing and deliver a total dose that would be efficacious. The scientific hypothesis is that the same total dose could be given safely with tolerated declines in Hb over a longer duration by starting with a lower initial dose which is increased gradually over time. The ascending dose regimen would allow for a steady adjustment of the age distribution of RBCs by both slow primaquineinduced haemolysis and the resulting increased erythropoiesis. These results only concern ascending dose regimens given over 20 days. There are two reasons for this; first, adherence to long course regimens is likely to be poor, and second, the first relapses emerge from the liver about 14 days after starting treatment so the primaquine regimen has to provide sufficient drug to prevent the emergence or eliminate these parasites.
Definition
For practical purposes, an acceptable ascending dose regimen is defined as a monotonic increasing dose regimen satisfying the following conditions: (i) the total dose is >380 mg , (current tablet sizes do not allow for a regimen to provide 420 mg easily to all adult patients—and 380 mg is considered to give similar efficacy); (ii) every increment is a multiple of 2.5 mg; (iii) the minimum adult daily dose is 5 mg; and (iv) the maximum adult daily dose is 30 mg.
The optimal ascendingdose regimen is defined as the one resulting in the slowest haemolysis, where the rate of haemolysis is penalized by the squared gradient. The optimisation problem is nonconvex for all ascending dose regimens, so the solution is approximated using a greedy search algorithm. An estimated optimal dosing regimen satisfying the criteria defined above is shown in Figure 6, plotted in red (left: haemolytic effect; right: daily dosing of the ascending regimen). This was found using the median Bayesian posterior parameter estimates and a doseresponse relationship taken from a linear interpolation of all points in Figure 5 (left plot). In blue is a simplified version of this ascendingdose regimen, broken into four 5day cycles at a fixed dose. The resulting haemolysis from the blue regimen is greater, and the drops in haemoglobin are more irregular (left plot). Video 1 in the supplementary materials illustrates the red blood cell dynamics over the course of this regimen.
Although intuitively one might think that starting with a lower dose was safer, such a regimen can be in fact worse. This is shown in Figure 5 which compares the haemolysis resulting from four regimens given in Table 1. Regimen D delivers too little primaquine at the start (observe very small decreases in haemoglobin concentration) with a reticulocyte response that is too weak to render the RBCs 'resistant' to primaquine; the necessity to increase the PQ dose too fast to compensate for the slow start and to deliver an efficacious total dose results in a large drop in Hb on day 22.
Discussion
Primaquine is widely recommended for the radical cure of vivax malaria but it is often not given because testing for G6PD deficiency is not widely available outside large centres. This has deleterious consequences for vivax malaria affected communities because it is the multiple relapses of vivax malaria from liver hypnozoites that cause substantial morbidity.
Seminal research conducted over 50 years ago characterized the biology of oxidant haemolysis caused by primaquine and provided an alternative once weekly regimen for patients who were G6PDd based on controlled haemolysis. This was shown to be safer in adult subjects with the 'mild' African ${A}^{}$ variant of G6PDd, but was recommended for all G6PDd variants with variable adoption by countries since. In some countries (for example, Iran) it is the standard radical treatment for all patients. The safety and effectiveness of the high dose weekly regimen have been studied little over the past five decades.
Uncomplicated malaria treatment recommendations are usually a tradeoff between dosing precision and operational feasibility. A regimen which is long or complicated may be adhered to poorly. In this particular case it must also be able to prevent or suppress relapsing P. vivax or P. ovale parasites which begin to emerge from the liver as early as two weeks (becoming patent about one week later) in SE Asia and Oceania. This modelling exercise, based on all available data, sought to devise a primaquine regimen which would be safer in G6PDd patients, and, therefore, might be deployed without G6PD testing. It was calibrated against recent data in Cambodian patients most of whom had the Viangchan G6PDd variant. Thus, the model predictions of the degree of haemolysis and the tolerability and safety profile would be expected to hold for variants with similar or less severe enzyme abnormalities, but it would not necessarily hold for more severe variants such as G6PDd Mediterranean where more clinical research is required.
Under all circumstances, the ascending regimen proposed here would be expected to be safer than the current 14 day regimens in G6PDd hemizygous males and homozygous females, especially the 0.5 mg/kg regimen needed for frequent relapsing P. vivax. This is clinically relevant also for female heterozygotes. Even with current rapid testing methods (for example, fluorescent spot test and RDTs) which generally detect patients with $\le 30$% normal G6PD activity, the haemolytic risk in heterozygote females, who may be classified erroneously as 'G6PD normal', could still be substantial. Up to $\approx 70$% of their erythrocytes may be G6PD deficient, and clinically significant haemolysis may result from daily higher dose primaquine regimens given to female heterozygotes (Chu et al., 2017).
Although this compartmental model of RBC dynamics is highly simplified, it reproduces the essential dynamics of the body’s response to primaquineinduced haemolysis in both healthy individuals and malaria patients. It can therefore help to guide the design of a Phase I study to evaluate its predictions, and thereby develop an optimal ascending dose regimen of PQ. An adaptive design protocol has been developed to test the simplified regimen (A) in G6PDd Mahidol healthy volunteers. A study in healthy G6PDd volunteers is essential to characterise the haemolytic response. Data from such a study can then be used to determine an optimal regimen which would then be tested for safety, and efficacy (that is, radical cure) in vivax malaria patients in a Phase II (that is, to define the PKPD relationship in patients). Whether patients would adhere sufficiently to a longer regimen is an important operational concern so the optimised regimen would then need to be assessed for safety and effectiveness in larger field trials.
This use of mathematical modelling such as this could also be readily applied to the slowly eliminated 8aminoquinoline tafenoquine, currently being tested for safety and efficacy in humans (Beck et al., 2016). Tafenoquine has the great advantage of being administered as a single dose for radical cure due to its long terminal elimination halflife. However, this means it could be dangerous in G6PD deficiency. Whereas the rapidly eliminated primaquine can be stopped if there is significant haemolysis, limiting the haemolytic effect, the haemolytic effect of the slowly eliminated tafenoquine cannot be readily reversed and so haemolysis will continue until all susceptible red cells are destroyed. Combined regimens for G6PDd patients in which primaquine is given initially to induce controlled haemolysis followed by tafenoquine might be possible, and would allow shorter total treatment durations. Data on the Hb response to different doses of tafenoquine would be necessary to calibrate the model.
The results of this study show how care will need to be taken when designing an ascending primaquine dose regimen in order to minimize falls in haemoglobin. This is shown by the toxic regimen D (Table 1 and Figure 7). This gives an insight into the ‘memory’ property of the ascending dose regimens. The effect of a 30 mg dose will entirely depend on which doses were given during the previous days. In conclusion, these results suggest that an ascending PQ dosing regimen for vivax radical cure might be well tolerated and effective in mild or moderately severe G6PDd variants. These predictions should now be tested in an adaptive phase I study.
Materials and methods
Mathematical model
The structure of the model of red cell dynamics is similar to the compartmental model developed by (Savill et al., 2009). RBC dynamics are simulated by tracking the age distribution of the red blood cells in hourly blocks. The homeostatic dynamics, which maintain the number of red blood cells or haematocrit at a steady state are straightforward. At steady state, approximately 0.83% of RBCs are replaced each day and 1% of RBCs in the circulation are reticulocytes. Severe acute anaemia has two consequences in the bone marrow. Reticulocytes are released into the circulating blood at an earlier age and with increased erythropoiesis normoblasts may be released into the circulation (reported as nucleated RBCs) (Hillman, 1969). Previous iron turnover studies in humans following phlebotomy suggest sigmoid relationships for both of these processes (Hillman, 1969).
Compartmental model of RBC dynamics
Request a detailed protocolThe steady state haemoglobin concentration is denoted $\mathrm{H}\mathrm{b}}^{\ast$. The time at which reticulocytes are released into circulation is a function of the haemoglobin concentration at time t, denoted Release ($\mathrm{H}\mathrm{b}}_{\mathrm{t}$). At steady state, when $\mathrm{H}\mathrm{b}}_{\mathrm{t}}={\mathrm{H}\mathrm{b}}^{\ast$, the reticulocytes mature in the bone marrow for 3.5 days (that is, Release $\{{\mathrm{H}\mathrm{b}}^{\ast}\}=3.5$) and then spend approximately one day in the circulation before becoming erythrocytes. In anaemia, release can occur after only one day (then the cells are reticulocytes for 3.5 days in the circulation). These relationships were shown from plasma iron turnover studies following phlebotomy (Hillman, 1969).
By modelling the number of circulating RBCs directly, it is possible to compute the haemoglobin concentration at each time point in the simulations. The steady state number of RBCs corresponds to steady state haemoglobin and steady state haematocrit. The relationship between haematocrit and haemoglobin is assumed to be linear (Lee et al., 2008). At steady state the body produces approximately ${10}^{8}$ RBCs per hour. In the model, this is used as the baseline production quantity, represented by a production factor of $\rho =1$. In extreme anaemia this can be increased fivefold or more, for example, $\rho \ge 5$. Both Release ($\mathrm{H}\mathrm{b}}_{\mathrm{t}$) and $\rho ({\mathrm{H}\mathrm{b}}_{t})$ are modelled as sigmoid functions:
with $\lambda ({\mathrm{H}\mathrm{b}}^{\ast})$ given by:
where $\mathrm{H}\mathrm{b}}_{t$ is the haemoglobin concentration at time $t$ (or equivalently the haematocrit); ${\rho}^{\mathrm{max}}$ is the maximum fold increase in steady state RBC production; $\mathrm{H}\mathrm{b}}^{\ast$ is the steady state haemoglobin concentration; $\mathrm{H}\mathrm{b}}_{50}^{\rho$ is the midpoint of the $\rho $ sigmoid function (the haemoglobin concentration at which production is equal to ${\rho}^{\mathrm{m}\mathrm{a}\mathrm{x}}/2;$$k$ the Hill coefficient which regulates the steepness in response to perturbations in haemoglobin levels; $\mathrm{H}\mathrm{b}}_{50}^{\mathrm{R}$ the midpoint of the release function.
For simplicity, it is assumed that in the normal healthy state all red blood cells live exactly 120 days as erythrocytes, and therefore these two functions are sufficient to model the feedback loops which regulate perturbations to haemoglobin levels.
The following class of functions is used to model the red cell agedependent primaquineinduced haemolysis:
where ${T}_{\mathrm{min}}$ is the age of the youngest red blood cells lysed by primaquine. This parameter varies as a function of the degree of G6PD deficiency (determined by the genetic variant of G6PDd). The data from (Pannacciulli et al., 1965) suggest cells as young as 16 days can be lysed with a daily primaquine dose of 30 mg in the severe Mediterranean variant, whereas in the less severe African A variant haemolysis appears confined to cells older than 50 days (Beutler et al., 1954a). The steepness of this ‘killing function’ is regulated by the parameter $d$, with smaller values of $d$ giving a sharper drop in haemoglobin levels.
Depending on the degree of severity of G6PD deficiency, haemolysis will be observed more or less quickly after the first dose of primaquine. To simulate this effect, a time lag component is added to the model, the value of which will depend on the genetic variant of G6PDd. The time lag component, denoted ${T}_{lag}$, reduces the total haemolytic effect (as given in Equation 1 for the first ${T}_{lag}$ hours after the first dose of primaquine. A glossary of all parameters of the model alongside their units and interpretation is given in Table 2.
Available data on primaquine induced haemolysis are sparse (these data are reviewed in detail in the next section). Thus, many of the free parameters such as ${T}_{lag}$ and ${T}_{\mathrm{min}}$ are fixed using expert opinion and their impact on the modelling evaluated by a sensitivity analysis shown in appendix Modelling primaquineinduced haemolysis in G6PD deficiency. A primary goal of this analysis was to parameterize the relationship between primaquine dose and haemolysis for a given severity of G6PD deficiency. Thus the doseresponse curve varied with the different genetic mutations. A preliminary analysis (Figure 5) of the available historical data (Figure 2) suggests that the relationship between dose and agedependent haemolysis as parameterised by the function in Equation 1 is logarithmic:
where $x$ are individual covariates of importance such the G6PDd variant and sex (the deficiency is Xlinked) and $PM{Q}_{\mathrm{d}\mathrm{o}\mathrm{s}\mathrm{e}}$ is the dose of primaquine.
Inputs and outputs of model
Request a detailed protocolThe goal of this compartmental model was to simulate haemolysis following primaquine administration to a G6PDd individual over a fixed period of time. In order to do this, the model needs as inputs the following elements.
The age distribution of RBCs both in the bone marrow (normoblasts and reticulocytes) and in circulation (reticulocytes and erythrocytes). This is represented as a vector of counts of RBCs for each age group discretized into hourly blocks.
The number of hours for which to simulate the model forward in time. This includes a dosing schedule (binary vector discretized in hourly blocks where 1 represents drug in the body and 0 represents drug absent). In all the simulations shown, the drug schedules are designed in multiples of 24 hr blocks (the terminal elimination halflife of primaquine is $\approx 5$ hours). In this manner it is possible to simulate weekly dosing and daily dosing.
The dosing level of the drug (this is defined by the value of parameter $d$ from Equation 1).
The simulations output the haemoglobin concentration over time, the reticulocyte percentage over time, and the final age distribution of RBCs in both the bone marrow and in circulation.
The structure of the model also allows for the simulation of the effect of malaria within an individual by altering the age distribution as P. vivax invades young erythrocytes exclusively. It is assumed that healthy individuals with no history of haemolytic events in the last 4 months will have a uniform age distribution of RBCs. Simulations in P. vivax infections could be done by shifting the age distribution. Model code is available in Source code 1. Supplementary file 1 provides posterior MCMC samples from the model.
Data
Historical studies
Request a detailed protocolAlthough primaquine was first tested in humans in 1944 and approved by the US FDA in 1952, there are very few precise data on its haemolytic effect in G6PDd individuals. Most of the studies only present sparse data, often limited to summary statistics, and with small sample sizes ($n\approx 2$). This section presents an analysis of the most informationrich studies conducted over the last 60 years and shows how they can be used to design informative prior distributions for the compartmental model. An exhaustive review of all available data on the haemolytic effects of primaquine has been reported recently by Recht et al. (2014). There are however only five studies which present useful data on falls in haemoglobin concentrations over time. Throughout this section, for consistency in Figures 1–4, data are shown in green and the model fits/predictions in red.
The data extracted from past studies in healthy volunteers are shown in Figures 2,3. These figures compare least squares model fits with the data. The studies span four different variants of G6PD deficiency: Mediterranean, Mahidol, Viangchan (The Khmer soldiers whose data are shown in Figure 2 are incorrectly referred to by Everett et al. (1977) as G6PDd Mahidol. They are in fact most likely to be G6PDd Viangchan, as the variant was discovered some ten years later [1988]) and African ${A}^{}$. These can be categorised generally as severe, moderate, moderate and mild variants of G6PDd, respectively.
The deterministic compartmental model reproduces the essential patterns of the observed dynamics, namely, a more rapid decrease in haematocrit at the start of the regimen which then slows as the haemolysis becomes selflimiting. However, some aspects cannot be reproduced by the model. The bottom left plot of Figure 2 shows a sharp drop in haematocrit after the second round of primaquine administration. This is surprising, although the fact that it is outside the administration period may be an error in plotting in the original paper. The compartmental model cannot reproduce such a marked second drop after an initial fall of over 50% from baseline. In Figure 3, for the bottom two plots, the model predicts a faster initial decrease than observed for a similar nadir drop.
These data were used to select a plausible range of values for the parameter $d$ in Equation 1. This parameter governs the agedependent haemolysis at different dosing levels.
Weekly highdose primaquine in G6PDd viangchan
Request a detailed protocolThe largest and most recent study of primaquineinduced haemolysis in G6PD deficiency is from Kheng et al. (2015). In this study, 75 Cambodian patients with vivax malaria were given primaquine 0.75 base mg/kg weekly for 8 weeks. Of these, 17 were G6PDd Viangchan (14 homozygous males and three heterozygous females) and 1 was a homozygous male with G6PDd Canton. Haemoglobin concentrations and reticulocyte counts were measured on days 0,1,2,3 and subsequently each week, before the next dose of primaquine. One patient had a marked drop of haemoglobin falling to 7.5 g/dL and required a blood transfusion. The data are shown in Figure 4. The variation in the measurement of haemoglobin was approximately 1 g/dL. For reticulocyte counts, the measurement error is much greater when done by microscopy (that is, counting per 100 red blood cells).
Model fitting using MCMC
Request a detailed protocolThe historical (Figures 2,3) and Kheng (Figure 4) data were used to fit the compartmental model; the former were used to select suitable prior distributions for parameters. Bayesian model fitting via MCMC was then applied to the data from the weekly highdose primaquine in Cambodia (Kheng et al., 2015). The likelihood of the parameters is defined by a deterministic simulation from the compartmental model for a given dosing regimen and assumes both the haemoglobin levels and reticulocyte counts are observed with Laplace distributed errors (assuming Laplace errors is equivalent to minimizing absolute deviation). A Bayesian hierarchical structure was used for the steady state haemoglobin and the maximum increase in the production of red blood cells. This makes the assumption that each patient in the study is characterized by an individual steady state haemoglobin concentration $\mathrm{H}\mathrm{b}}^{\ast$ and a maximum production capacity ${\rho}^{\mathrm{max}}$ drawn from a population distribution (normal distributions in both cases). Weekly informative priors were used for all parameters and the posterior distribution was estimated using MCMC with a MetropolisHastings proposal. Details of prior distributions and histograms of posterior distributions, together with convergence diagnostics and summary statistics are in the Appendix 1, ‘Structure of hierarchical model and MCMC diagnostics'.
Appendix 1
Structure of hierarchical model and MCMC diagnostics
The analysis from section ‘Model fitting’ assumes that the parameters $d,k,{\mathrm{H}\mathrm{b}}_{50}^{\rho},{\mathrm{H}\mathrm{b}}_{50}^{R}$ are defined on a population level (same for all individuals in the study). The parameters $\mathrm{H}\mathrm{b}}^{\ast},{\rho}^{max$ are defined on the individual level but drawn from a population distribution. Thus the ${i}^{th}$ patient is characterized by:
We set informative priors on the hyperparameters $\theta}_{\mathrm{H}\mathrm{b}},{\theta}_{\rho$, namely:
where ${N}_{trunc}$ is a truncated normal distribution with upper and lower values given by $\mathrm{min}$ and $\mathrm{max}$. The variance hyperparameters $\sigma}_{\mathrm{H}\mathrm{b}}^{2$ and ${\sigma}_{\rho}^{2}$ are given flat priors. Figure 2 shows the posterior distributions for all the parameters and hyperparameters of the model. The population level parameters are given the following informative priors:
Convergence of the Metropolis Hastings algorithm was done by running four independent chains and computing the GelmanRubin statistic.
As an independent check for model fit, Figure 3, compares the estimated individual steady state haematocrits with the data on individual mean corpuscular volume (MCV) on day 0. The relationship between these fitted estimates and the data from the study is approximately linear for all patients (G6PDd and G6PD normal). The MCV is a poor marker of anaemia; high values are associated with reticulocytosis, folate and B12 deficiency and low values with iron deficiency, and thalassaemia both of which are very common in tropical areas. A correlation with steady state haematocrit would suggest that the model is indeed estimating the correct quantities. Therefore if the model estimates of the steady state haematocrit are correct, this should in theory correlate well with the baseline MCV and suggest the model is converging to the correct quantities.
Sensitivity analysis
Figure 4 shows the effect of the individual parameters of the model (except the killing parameter $d$, the main parameter if interest in the model) on the haematocrit and reticulocyte response when all other parameters are held fixed.
To illustrate the individual effects of the parameters, we take as example a daily dosing regimen with identical dose (that is, same value of the parameter $d$) for 30 days. This is chosen to approximate the drop in haemoglobin observed in G6PDd African ${A}^{}$ studies (Figure 2). We fix $d=0.05$.
Some notable points:
The maxfold increase (${\rho}^{\mathrm{max}}$) of RBCs in the bone marrow has the most impact on the reticulocyte response after day 9. We assume that a twofold increase in production is the minimum viable response to anaemia (shown by the red lines). Hillman (1969) estimates the maxfold increase in healthy males as approximately 5fold.
The analysis of the data from Kheng et al. (2015) gives very similar values for the population distribution of $\mathrm{H}\mathrm{b}}_{50}^{\rho$ as that estimated by Hillman (1969) (see Figure 2 in Hillman [1969]).
References

Mitigation of the haemolytic effect of primaquine and enhancement of its action against exoerythrocytic forms of the chesson strain of piasmodium vivax by intermittent regimens of drug administration: a preliminary reportBulletin of the World Health Organization 22:621.

Estimation of the antirelapse efficacy of tafenoquine, using Plasmodium vivax GenotypingJournal of Infectious Diseases 213:794–799.https://doi.org/10.1093/infdis/jiv508

The hemolytic effect of primaquine. IV. the relationship of cell age to hemolysisThe Journal of Laboratory and Clinical Medicine 44:439–442.

The hemolytic effect of primaquine. III. A study of primaquinesensitive erythrocytesThe Journal of Laboratory and Clinical Medicine 44:177–184.

The hemolytic effect of primaquine. VII. biochemical studies of drugsensitive erythrocytesThe Journal of Laboratory and Clinical Medicine 45:286–295.

The hemolytic effect of primaquine and related compounds: a reviewBlood 14:103–139.

Validation of the rapid test carestart(tm) G6PD among malaria vivaxinfected subjects in the brazilian amazonRevista da Sociedade Brasileira de Medicina Tropical 49:446–455.https://doi.org/10.1590/0037868201342016

The haemolytic effect of a single dose of 45 mg of primaquine in G6PD deficient thaisJournal of the Medical Association of Thailand = Chotmaihet thangphaet 55:631–638.

The hemolytic effect of primaquine. II. the natural course of the hemolytic anemia and the mechanism of its selflimited characterThe Journal of Laboratory and Clinical Medicine 44:171–176.

Hemoglobin E and glucose6phosphate deficiency in the khmer air force (Cambodia)The American journal of tropical medicine and hygiene 26:597–601.

A long neglected world malaria map: plasmodium vivax endemicity in 2010PLoS Neglected Tropical Diseases 6:e1814.https://doi.org/10.1371/journal.pntd.0001814

Characteristics of marrow production and reticulocyte maturation in normal man in response to anemiaJournal of Clinical Investigation 48:443–453.https://doi.org/10.1172/JCI106001

G6PD deficiency: global distribution, genetic variants and primaquine therapyAdvances in Parasitology 81:133–201.https://doi.org/10.1016/B9780124078260.000047

The course of experimentally induced hemolytic anemia in A primaquinesensitive caucasian. a case studyBlood 25:92–95.

Vivax malaria: neglected and not benignThe American Journal of Tropical Medicine and Hygiene 77:79–87.

BookSafety of 8aminoquinoline antimalarial medicinesGeneva: World Health Organization.

Nature of hemolytic crises and the fate of G6PD deficient, drugdamaged erythrocytes in SardiniansNew England Journal of Medicine 276:1339–1344.https://doi.org/10.1056/NEJM196706152762402

Assessment of PointofCare diagnostics for G6PD deficiency in malaria endemic rural eastern indonesiaPLOS Neglected Tropical Diseases 10:e0004457.https://doi.org/10.1371/journal.pntd.0004457

Quantitative analysis of mechanisms that govern red blood cell age structure and dynamics during anaemiaPLoS Computational Biology 5:e1000416.https://doi.org/10.1371/journal.pcbi.1000416
Article and author information
Author details
Funding
Wellcome
 James Watson
 Walter RJ Taylor
 Nicholas J White
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
This study was part of the Wellcome Trust Mahidol University Oxford Tropical Medicine Research Programme funded by the Wellcome Trust. We thank Nick Savill for sharing his code. We thank Lisa White and Ben Cooper for help with designing the model.
Copyright
© 2017, Watson et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 1,903
 views

 329
 downloads

 43
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading

 Computational and Systems Biology
 Developmental Biology
Understanding the principles underlying the design of robust, yet flexible patterning systems is a key problem in developmental biology. In the Drosophila wing, Hedgehog (Hh) signaling determines patterning outputs using dynamical properties of the Hh gradient. In particular, the pattern of collier (col) is established by the steadystate Hh gradient, whereas the pattern of decapentaplegic (dpp), is established by a transient gradient of Hh known as the Hh overshoot. Here we use mathematical modeling to suggest that this dynamical interpretation of the Hh gradient results in specific robustness and precision properties. For instance, the location of the anterior border of col, which is subject to selfenhanced ligand degradation is more robustly specified than that of dpp to changes in morphogen dosage, and we provide experimental evidence of this prediction. However, the anterior border of dpp expression pattern, which is established by the overshoot gradient is much more precise to what would be expected by the steadystate gradient. Therefore, the dynamical interpretation of Hh signaling offers tradeoffs between

 Computational and Systems Biology
 Neuroscience
Animal behaviour alternates between stochastic exploration and goaldirected actions, which are generated by the underlying neural dynamics. Previously, we demonstrated that the compositional Restricted Boltzmann Machine (cRBM) can decompose wholebrain activity of larval zebrafish data at the neural level into a small number (∼100200) of assemblies that can account for the stochasticity of the neural activity (van der Plas et al., eLife, 2023). Here, we advance this representation by extending to a combined stochasticdynamical representation to account for both aspects using the recurrent temporal RBM (RTRBM) and transferlearning based on the cRBM estimate. We demonstrate that the functional advantage of the RTRBM is captured in the temporal weights on the hidden units, representing neural assemblies, for both simulated and experimental data. Our results show that the temporal expansion outperforms the stochasticonly cRBM in terms of generalization error and achieves a more accurate representation of the moments in time. Lastly, we demonstrate that we can identify the original timescale of assembly dynamics by estimating multiple RTRBMs at different temporal resolutions. Together, we propose that RTRBMs are a valuable tool for capturing the combined stochastic and timepredictive dynamics of largescale data sets.