Nonlinear mixed effects modeling of gametocyte carriage in patients with uncomplicated malaria
© Distiller et al; licensee BioMed Central Ltd. 2010
Received: 28 May 2009
Accepted: 26 February 2010
Published: 26 February 2010
Gametocytes are the sexual form of the malaria parasite and the main agents of transmission. While there are several factors that influence host infectivity, the density of gametocytes appears to be the best single measure that is related to the human host's infectivity to mosquitoes. Despite the obviously important role that gametocytes play in the transmission of malaria and spread of anti-malarial resistance, it is common to estimate gametocyte carriage indirectly based on asexual parasite measurements. The objective of this research was to directly model observed gametocyte densities over time, during the primary infection.
Of 447 patients enrolled in sulphadoxine-pyrimethamine therapeutic efficacy studies in South Africa and Mozambique, a subset of 103 patients who had no gametocytes pre-treatment and who had at least three non-zero gametocyte densities over the 42-day follow up period were included in this analysis.
A variety of different functions were examined. A modified version of the critical exponential function was selected for the final model given its robustness across different datasets and its flexibility in assuming a variety of different shapes. Age, site, initial asexual parasite density (logged to the base 10), and an empirical patient category were the co-variates that were found to improve the model.
A population nonlinear modeling approach seems promising and produced a flexible function whose estimates were stable across various different datasets. Surprisingly, dihydrofolate reductase and dihydropteroate synthetase mutation prevalence did not enter the model. This is probably related to a lack of power (quintuple mutations n = 12), and informative censoring; treatment failures were withdrawn from the study and given rescue treatment, usually prior to completion of follow up.
Gametocytes are the sexual form of the malaria parasite and the main agents of transmission. Malaria is diagnosed based on the presence of asexual parasites and most anti-malarials act on these stages with the aim of eliminating them. However, of vital importance to control (and eventually eliminate) malaria and prevent the spread of drug resistance is the prevention and/or curtailment of gametocytemia [1, 2].
While transmission of the malaria parasite is dependent on various factors such as vector characteristics, host susceptibility and climatic conditions, it is host infectivity (the probability of a mosquito becoming infected after biting an infected person) that plays a crucial role in the onward spread of an individual infection [3–5]. There are various factors that influence human host infectivity to mosquitoes and the density of gametocytes appears to be the best single measure available .
Despite the obvious importance of gametocyte carriage for both the transmission of malaria and the spread of resistance, little research has been conducted modeling gametocyte densities, and by extension host infectivity, directly. Rather, models generally rely on the more frequently measured asexual parasite densities from which gametocyte carriage is estimated by either making a simple assumption (such as a constant gametocyte switching rate), or by using a more advanced model such as using lagged values based on asexual parasite densities [6, 7].
Yet the relationship between gametocytes and asexual parasites is not straightforward. While one would expect a higher asexual parasite load to be associated with a higher sexual parasite load, there are also other factors that could be associated with low asexual densities such as immunity, treatment or a long duration of infection. As some of these factors would also increase gametocyte carriage, the relationship in question is often confounded in a clinical setting [2, 8].
Repeated measures: Patients are seen on up to nine occasions during a 42-day period and gametocyte densities are measured on seven of these occasions.
Study design: Studies are designed to monitor asexual parasite density as an indicator of therapeutic efficacy. Gametocyte emergence lags asexual parasite emergence and hence the intervals for measurement of gametocytemia are not ideal.
Informative censoring for treatment failures: Patients are withdrawn and given rescue treatment if parasites do not clear or if they re-emerge, resulting in censoring of gametocyte distributions. In this way the failure to clear asexual parasites is a competing event to the measuring of gametocytes, with the result that the gametocyte distributions are not seen or are incomplete for patients who are resistant to the anti-malarial being studied and hence are withdrawn due to treatment failure.
One way of dealing with the bias due to informative censoring is, or could be, by imputation of the missing parts of the gametocyte distributions. To do this, a method that directly models the gametocyte density-time relationship for different groups of patients is needed.
Data for this research came from open label clinical trials of patients treated with sulphadoxine-pyrimethamine (SP) for uncomplicated symptomatic falciparum malaria at the Naas and Mangweni clinics (Mpumalanga province, South Africa) and Namaacha and Bela Vista Clinics (Maputo province, Mozambique). The study design focused on measuring the clinical and asexual parasitological response to treatment as recommended by the World Health Organization. Gametocyte densities were thus measured on days 0, 3, 7, 14, 21, 28, and 42 following treatment.
Of the 447 patients enrolled in these studies [9–11], a subset of 103 patients who had no gametocytes pre-treatment and who had at least three non-zero gametocyte density measurements (including 56 patients with at least four non-zero readings) were included in this analysis. Patients with gametocytes pre-treatment (n = 23) were excluded so that the analysis could focus on the development of gametocytes during the primary infection, rather than pre-existing gametocytes potentially from previous infections. Gametocyte density was counted on thick blood smears against 1,000 leukocytes assuming 8,000 leukocytes per microlitre (μL). This explains why the lowest detectable density was 8 per μL as the smallest possible count would have been one gametocyte per 1,000 leukocytes.
Patients who failed treatment (in that they failed to clear their asexual parasites) were withdrawn from the study to be given rescue treatment and effectively were lost to follow up. Due to the time lag between asexual parasites and sexual gametocytes, the distribution of gametocyte densities for these patients were often truncated or, in some cases, the patients were withdrawn prior to the emergence of gametocytes. This resulted in informative censoring of gametocyte density distributions for treatment failures.
The original gametocyte readings were logged to the base two for the modeling. It is common to log a measurement that is characterized by long right tails and where it is more meaningful to express increases in a multiplicative rather than additive manner. Using a logarithmic scale with base two allowed a one unit change in the response to be viewed as a doubling or halving. For example, the gametocyte elimination half-life is the time it takes for the density of gametocytes to halve.
The main problem with the data was the high prevalence of zero measurements and their associated uncertainty. Attempting to include many zeros in the data when trying to solve nonlinear mixed effect models led to estimation difficulties. Therefore there were two broad treatments of these zeros. In both cases the zero reading at day 0 was left at zero and patients exhibiting gametocytes on day 0 were excluded. In the first dataset (hereafter referred to as Data 1) all subsequent zeros were put to missing, while in the second instance (Data 2) they were changed to the detectable limit of eight gametocytes per microlitre (μL). Histograms were examined showing the overall frequency distributions and found that the shape of the distribution does not change qualitatively when we change the eight to a four or a twelve.
These studies were approved by the Research and Ethics Committees of the University of Cape Town, the Mpumalanga Department of Health and the Mozambican Ministry of Health prior to commencement. These clinical studies were conducted in accordance with the South African Clinical Trials Guidelines  and the principles laid down by the World Health Assembly of 1975 on Ethics in Human Experimentation and the Helsinki Declaration. Written Informed consent was obtained from each subject or their guardian. In the case of illiteracy, an "X" mark was accepted as documentation of consent in the presence of an independent literate witness who also signed the consent form. Study subject identities were kept confidential. Those who chose not to participate in this study were treated with sulphadoxine-pyrimethamine but without routine follow up. At the time of the study, sulphadoxine-pyrimethamine and chloroquine monotherapies were the malaria treatment policies in Mpumalanga Province, South Africa and Mozambique, respectively.
Empirical patient categories
group 1: this was the largest patient category and included patients that experienced an initial increase followed by a decrease that did not increase again.
group 2: the second group was for patients who experienced a second increase in gametocytes after the initial rise had started going down but the increase was slight i.e. less than a twofold increase.
group 3: the third category included those patients that experienced a marked second increase in gametocytes after the initial rise had started going down. The marked increase was defined as a change of at least twofold.
Univariate analysis of determinants of patient category.
1 (n = 74)
2 (n = 14)
3 (n = 15)
Parasite Clearance Time
Logged Parasite Density
The gametocyte density-time profiles were typically nonlinear and were further characterized by both intra-individual variability due to repeated measures and inter-individual variability, which made nonlinear mixed effects modeling the appropriate statistical methodology. The modeling was conducted using R 2.2.0  and in particular the nlme package .
Finding an underlying (nonlinear) function with a common structure that could be generating the data: y ij = f (t ij , β i ) + e ij where y ij is the observed response for the ith subject at the jth time point, f (t ij , β i ) refers to a function f whose structure is common to all patients with t ij representing time and parameters β i representing the value for the function's parameter(s) for the ith person, and e ij refers to the error or residual for the ith person at the jth time point and e ij ~ N (0, V i ) with V i being the intra-individual covariance matrix.
Inter-individual variability is modeled by adding random effects to mean parameters to allow individual specific responses, β i = β + b i where β i is subject i's parameter value for the function and can be thought of as some random deviation (b i also know as a random effect) from the mean parameter value (β). It is assumed that b i ~ N (0, D) where D is a (k × k) dispersion matrix for the k random effects, and that the b i are independent of each other and the e ij .
Co-variates can also be used to explain systematic variability in the actual parameter values between subjects: β i = Aβ + b i where the matrix A is a design matrix formed from one or more co-variates.
With regards to the random effects, typically emphasis is placed not on the actual value of the random effects (b i 's) but rather the size of their variability. The covariance matrix of the random effects is specified with the model and includes both diagonal elements (variances of the random effects for particular parameters) and off-diagonal elements (covariance of the random effects between different parameters). More complicated structures are obviously more difficult to solve.
Various subsets of the data were modeled (not all shown here) and these models were compared by looking at both the Akaike Information Criteria's (AIC's), and at how robust the functions were across different datasets. The AIC is the appropriate statistic for comparing different models fitted by maximum likelihood. The aim was to find a function that fitted the data well, that made biological sense, and that was not overly sensitive to either the treatment of the zeros or the inclusion of different patient categories.
Finding the function
In total six different nonlinear functions were initially fitted to the different datasets. These functions included classic pharmacokinetic functions such as the one and two compartment models, the critical exponential (a variation of the double exponential), a reduced form of the critical exponential (hereafter called the "modified critical exponential"), and wave-like functions, such as the fourier and double fourier. Models based on different nonlinear functions were compared using the AIC statistic. Mixed effect modeling started with the most general and complex variance-covariance structure i.e. random effects for all function parameters were included and nonzero correlations amongst all of them were allowed. If this complex structure could be estimated, the results were examined for signs of possible improvement/simplification, i.e. non-significant and/or unnecessary random effects were dropped. If the model with the most general variance-covariance structure could not be estimated, a simpler structure was attempted. The resulting nested models were compared using likelihood ratio (LR) χ2 statistics.
Use of the LR test in this strategy meant that Maximum Likelihood (ML) was used as opposed to Restricted Maximum Likelihood (REML) in the estimation process; any change to either the fixed or the random effects structure affected the REML likelihood and hence different structures could not be compared in this way if REML was used.
This phase examined if any co-variates could improve the model and explain some of the individual variation in the model parameters, thereby reducing the random effects. As suggested by Pinheiro and Bates , the strategy was to fit a model without any co-variates and then plot the random effects from this model against the potential co-variates. Any systematic pattern in the relationship between a potential co-variate and a particular random effect suggested that the particular co-variate be included in the model. Co-variates were added to the model one at a time starting with the one with the most obvious pattern. The new model was fitted and the process repeated with the new set of random effects. In addition to the LR test for nested models, statistical significance of any fixed effect associated with a co-variate was examined and various interactions were also explored.
The prospectively defined co-variates used in this analysis included: site (Mpumulanga vs Maputo province as the reference), age in years, gender (males versus females), day 0 parasite density (logged to the base 10), mutations associated with SP resistance (quintuple versus 0-4 mutations) , treatment outcome (failure versus success), parasite clearance time (either 1, 2, 3 or 7 days), and the empirical patient category described above.
Results and Discussion
It was apparent that the AIC from the one compartment model was better than the other two functions for Data 1, but unfortunately the one compartment model did not perform well with Data 2. In Data 1 the zeros were dropped thereby allowing a smoother function to be fitted as the fitted curve did not need to reach down to the extremely small observations later on. When the zeros were replaced by a small number (in this case 8) in Data 2, the function then attempted to take these dips into account.
Critical exponential model & modified version.
Fixed Effects (std error)
Fixed at 0
Fixed at 0
Fixed at 0
Fixed at 0
For the critical exponential function fitted to Data 1, the estimated A and B parameters cancelled each other out and the lack of a zero or near to zero reading at the last time point resulted in a biologically implausible asymptotic level of approximately 5 (or 25 = 32 μL). For Data 2, both the A and B parameters were estimated to be essentially (or close to) zero.
The modified critical exponential function (where the A and B parameters were fixed to zero) was the final function chosen as most appropriate for this data. This decision was based on the clear stability that this structure demonstrated across various datasets (including a subset of patients that had at least four non-zero gametocyte density readings), particularly with regards to the parameter R. It is worth noting that the chosen function has only two fixed effect parameters and that the data did not seem able to support more complicated functions. For Data 2, fixing the A and B parameters to equal zero enabled us to estimate random effect components for the R parameter, and resulted in a better fitting model (AIC = 2776 vs 2877).
The modified critical exponential function is similar to an ordinary exponential decay curve . This is apparent when it is rewritten as y = (C * t) × exp(logRt)where logR can be substituted for -β1, thereby constraining R to fall between zero and one in order to give a negative value. The difference is that the first term (C * t) is not constant like β0 but rather depends on time. The ordinary exponential decay curve has the property that the ratio of the response at successive time points ( ) is equal across the range of the function. This property does not hold for the modified critical exponential function.
The second stage of modeling resulted in the co-variates age, site (Mpumulanga vs Maputo province), logged baseline parasite density, and empirical patient category (pcat) being included in the model. The final model formulation was thus as follows:
and: e ij ~ N (0, V i )
The Final Model.
Intercept (βC 0)
(1.8851 : 2.3303)
Mpm vs Map (βC 1)
(-0.7503 : -0.2472)
Pcat2 vs 1 (βC 2)
(-0.3810 : 0.2630)
Pcat3 vs 1 (βC 3)
(-0.7989 : -0.1832)
Intercept (βR 0)
(0.8685 : 0.9007)
Mpm vs Map (βR 1)
(0.0093 : 0.0213)
Age (years) (βR 2)
(0.0001 : 0.0003)
log 10 (pdens) (βR 3)
(0.0027 : 0.0095)
Pcat2 vs 1 (βR 4)
(0.0032 : 0.0179)
Pcat3 vs 1 (βR 5)
(0.0058 : 0.0210)
(0.3964 : 0.5944)
(0.0078 : 0.0130)
(-0.8479 : -0.5872)
(1.2864 : 1.4597)
Patients from Mpumulanga tended to reach lower gametocyte densities compared to patients from Maputo and had slower rates of clearance (as indicated by a negative value for βC 1and a positive value for βR 1, while older people and people with high pre-treatment asexual parasite densities tended to have higher gametocyte densities (Table 4 shows that the R parameter varied according to age and baseline parasite density level).
The site effect could reflect the level of immunity as malaria transmission in Maputo, Mozambique is much more intense (> 20-fold higher prevalence) than in Mpumulanga, South Africa - and delayed treatment seeking in partially immune patients may explain the faster increase as well as the faster clearance rate seen in Mozambique. The site effect could also reflect other unmeasured characteristics such as genetic differences.
Since 1900 it has been recognized that the frequency and density of asexual parasites decreases with increasing age in areas with stable malaria transmission as a result of partial immunity being acquired [6, 18]. Unfortunately there is no direct measure of this immunity, so age is often used as a proxy to measure the effect of acquired immunity on various clinical and parasitological outcomes (asexual parasite density, failure rates, severe malaria rates) [6, 19]. As immunity increases with age in areas of stable malaria transmission, lower immunity in children would result in children seeking treatment earlier than adults, explaining the slightly higher gametocyte densities in adults. As expected, higher asexual parasite densities at baseline were related to higher gametocyte densities.
Lastly, the model predicted that patients in category two, on average, had a higher and later peak than category one but with a very similar rate of increase and decrease. On the other hand, the model predicted that patients in the third category had a lower peak that occurred slightly later and then cleared more slowly, i.e. these patients would have carried gametocytes for longer. This made sense when one considered that patients in the third category had a second (and sometimes third) wave of gametocytes, and that the function used in the model coped with this by smoothing out these peaks, thereby resulting in a lower average peak and a slower clearance rate.
It was interesting that there was very little difference in the effect of the other co-variate parameters when excluding patient category. This suggested that the patient categories could be viewed as a proxy for characteristics that were not measured (ie immunity, duration of primary infection and degree of synchronicity), yet were also independent from the available co-variate information.
Host infectivity is a critical component in the transmission of malaria and the spread of anti-malarial drug resistance. While various factors are thought to influence host infectivity, the best indicator appears to be gametocyte density. The huge obstacle that drug resistance poses for the control of malaria together with the role that gametocytes play in the spread of such resistance makes it critical to examine gametocyte carriage more closely.
However, most malaria transmission models are based on asexual parasite distributions that infer gametocyte densities through assumed switching rates. Yet the relationship between asexual and sexual parasites is complex and so models that can cope with modeling gametocyte data directly are desirable. The model developed in this research aimed to make a contribution by modeling gametocytes directly. The study design from which this data was obtained (which is used for most clinical trials) used measurement points that were based around the expected clinical and asexual parasitological response to malaria treatment, and hence patients were seen fairly intensively during the first week of their infection and then less intensively later on, i.e. a gap of 14 days between the last two time points. Since the distribution of gametocytes lags that of asexual parasites the current design was not optimal for measuring gametocytes and consequently non-zero gametocyte densities were relatively sparse. This problem was exacerbated by the withdrawal of patients who failed treatment so that they could be given rescue treatment.
From a public health perspective, accurately describing gametocyte carriage (and hence infectivity) is of greater relevance at the population level than for individual patients and hence a population modeling approach was used to overcome the limitations of sparse data. The model developed with the modified critical exponential function seems able to encompass a wide range of variability in the curves, and appeared to be stable across the different datasets. The main use of the derived model would be to simulate curves that could be used to impute truncated/censored gametocyte distributions possibly using a multiple imputation method http://www.missingdata.org.uk.
The spread of resistant infections results from their increased gametocytaemia . Thus it was surprising that the resistant (quintuple) mutation variables did not enter the model. This may be related to power as there were very few patients (n = 12) included with highly resistant infections. Including those patients with gametoytes pre-treatment (many of which would be recrudescent infections) would improve power, as these patients had a greater tendency to carry the quintuple mutation [8/23 (35%) vs 12/103 (12%), pvalue = 0.001], though as stated above the objective was to model gametocyte densities during the primary infection. Furthermore, most treatment failures were rescued and withdrawn from the study prior to completion of follow up. This informative censoring makes it very difficult to accurately estimate these gametocyte density-time profiles in resistant infections. Either a questionable assumption (such as equating their curves to those cured) has to be made, or a different study design that monitors gametocytes after patients have been withdrawn is needed. The issue of informative censoring in gametocyte data is an avenue that needs to be explored further. More data containing patients with resistant infections would increase the power to detect a "mutation effect" and hence allow the estimation of curves for treatment failures.
Lastly, this research was limited to carriers of gametocytes. In order to model the dynamics of malaria transmission and the spread of resistance in a population, the determinants of gametocyte prevalence should be assessed as reported previously using zero-inflated models used by Barnes et al[16, 20].
The authors wish to acknowledge the financial support provided by the UNDP/World Bank/WHO Special Programme for Research and Training in Tropical Diseases (TDR) to the South East African Combination Antimalarial Therapy (SEACAT) Evaluation, under which this study was undertaken. We also thank the staff of the Maputo Provincial Department of Health in Mozambique and the Mpumalanga Malaria Control Programme in South Africa, and the staff at the study health facilities for their contributions to the data collected for this study.
- Targett G, Drakeley C, Jawara M, von Seidlein L, Coleman R, Deen J, Pinder M, Doherty T, Sutherland C, Walraven G, Milligan P: Artesunate reduces but does not prevent posttreatment transmission of Plasmodium falciparum to Anopheles gambiae. J Infect Dis. 2001, 183: 1254-1259. 10.1086/319689.View ArticlePubMedGoogle Scholar
- Barnes K, White NJ: Population biology and antimalarial resistance: The transmission of antimalarial drug resistance in Plasmodium falciparum. Acta Trop. 2005, 94: 230-240.View ArticlePubMedGoogle Scholar
- Draper CC: Observations on the infectiousness of gametocytes in hyper-endemic malaria. Trans R Soc Trop Med Hyg. 1953, 47: 160-165. 10.1016/0035-9203(53)90072-8.View ArticlePubMedGoogle Scholar
- Diebner HH, Eichner M, Molineaux L, Collins WE, Jeffery GM, Dietz K: Modelling the transition of asexual blood stages of Plasmodium falciparum to gametocytes. J Theor Biol. 2000, 202: 113-127. 10.1006/jtbi.1999.1041.View ArticlePubMedGoogle Scholar
- Killeen G, Ross A, Smith T: Infectiousness of malaria-endemic human populations to vectors. Am J Trop Med Hyg. 2006, 75: 38-45.PubMedGoogle Scholar
- Pongtavornpinyo W: Mathematical modelling of antimalarial drug resistance. PhD thesis. 2006, School of Tropical Medicine, University of LiverpoolGoogle Scholar
- Ross A, Killeen G, Smith T: Relationships between host infectivity to mosquitoes and asexual parasite density in Plasmodium falciparum. Am J Trop Med Hyg. 2006, 75 (2 Supplement): 32-37.PubMedGoogle Scholar
- Price R, Nosten F, Simpson J, Luxemburger C, Phaipun L, ter Kuile F, Vugt MV, Chongsuphajaisiddhi T, White N: Risk factors for Gametocyte carriage in uncomplicated Falciparum malaria. Am J Trop Med Hyg. 1999, 60: 1019-1023.PubMedGoogle Scholar
- Barnes K, Little F, Smith P, Evans A, Watkins M, White N: Sulfadoxine-pyrimethamine pharmacokinetics in malaria: Pediatric dosing implications. Clin Pharmacol Ther. 2006, 80: 582-596. 10.1016/j.clpt.2006.08.016.View ArticlePubMedGoogle Scholar
- Mabuza A, Govere J, Durrheim D, Mngomezulu N, Bredenkamp B, Barnes K, Sharp B: Therapeutic efficacy of sulfadoxine-pyrimethamine for Plasmodium falciparum malaria. S Afr Med J. 2005, 95: 346-9.PubMedGoogle Scholar
- Mabuza A, Govere J, Durrheim D, Mngomezulu N, Bredenkamp B, Barnes K, Sharp B: Therapeutic efficacy of sulfadoxine-pyrimethamine in uncomplicated Plasmodium falciparum malaria 3 years after introduction in Mpumalanga. S Afr Med J. 2001, 91: 975-8.PubMedGoogle Scholar
- Guidelines for good practice in the conduct of clinical trails in South Africa. Tech. rep., Department of Health, South Africa. 2000Google Scholar
- R Development Core Team: R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. 2005, [http://www.R-project.org]Google Scholar
- Pinheiro J, Bates D, DebRoy S, Sarkar D: nlme: Linear and nonlinear mixed effects models. 2006, [R package version 3.1-68.1]Google Scholar
- Pinheiro J, Bates D: Mixed-effects Models in S and S-PLUS. 2002, Statistics and Computing, New York: Springer-VerlagGoogle Scholar
- Barnes K, Little F, Mabuza A, Mngomezulu N, Govere J, Durrheim D, Roper C, Watkins B, White N: Increased Gametocytemia after Treatment: An Early Parasitological Indicator of Emerging Sulfadoxine-Pyrimethamine Resistance in Falciparum Malaria. J Infect Dis. 2008, 197: 1605-1613. 10.1086/587645.View ArticlePubMedGoogle Scholar
- Eastwood DC, Mead A, Sergeant MJ, Burton KS: Statistical modelling of transcript profiles of differentially regulated genes. BMC Mol Biol. 2008, 9: 10.1186/1471-2199-9-66.Google Scholar
- Baird JK: Age-dependent characteristics of protection v. susceptibility to Plasmodium falciparum. Ann Trop Med Parasitol. 1998, 92: 367-390. 10.1080/00034989859366.View ArticlePubMedGoogle Scholar
- Aron JL, May RM: The poulation dynamics of Malaria. Population Dynamics of Infectious Diseases, Population and Community Biology. Edited by: Anderson RM. 1982, London: Chapman and Hall, 139-179.View ArticleGoogle Scholar
- Barnes K, Durrheim D, Little F, Jackson A, Mehta U, Allen E, Dlamini S, Tsoka J, Bredenkamp B, Mthembu D, White N, Sharp B: Effect of artemether-lumefantrine policy and improved vector control on malaria burden in KwaZulu-Natal, South Africa. PLoS Med. 2005, 2: e330-10.1371/journal.pmed.0020330.PubMed CentralView ArticlePubMedGoogle Scholar
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.