ON THE ESTIMATION OF SEQUESTERED INFECTED ERYTHROCYTES IN PLASMODIUM FALCIPARUM MALARIA PATIENTS

The aim of this paper is to give a method for the estimation of total parasite burden of the patient and the rate of infection in a malaria’s intra-host model by using control theory tools. More precisely, we use an auxiliary system, called observer or estimator, whose solutions tend exponentially to those of the original model. This observer uses only the available measurable data, namely, the values of peripheral infected erythrocytes. It provides estimates of the sequestered infected erythrocytes, that cannot be measured by clinical methods. Therefore this method allows to estimate the total parasite burden within a malaria patient. Moreover, our constructed observer does not use the uncertain infection rate parameter β. In fact, we derive a simple method to estimate this parameter β. We apply this estimation method using real data that have been collected when malaria was used as therapy for neurosyphilis by the US Public Health Service.


(Communicated by Zhilan Feng)
Abstract.The aim of this paper is to give a method for the estimation of total parasite burden of the patient and the rate of infection in a malaria's intra-host model by using control theory tools.More precisely, we use an auxiliary system, called observer or estimator, whose solutions tend exponentially to those of the original model.This observer uses only the available measurable data, namely, the values of peripheral infected erythrocytes.It provides estimates of the sequestered infected erythrocytes, that cannot be measured by clinical methods.Therefore this method allows to estimate the total parasite burden within a malaria patient.Moreover, our constructed observer does not use the uncertain infection rate parameter β.In fact, we derive a simple method to estimate this parameter β.We apply this estimation method using real data that have been collected when malaria was used as therapy for neurosyphilis by the US Public Health Service.
1. Introduction.Malaria is a disease that causes at least one million deaths around the world each year, with ninety per cent among African children. the most dangerous type of malaria is caused by the most virulent species of the Plasmodium parasite: Plasmodium falciparum.Intra-host models describe the dynamics of the interaction of parasites (viruses, bacteria, protozoans, ...) within the host.These models have been used to understand the population dynamics and evolution of Plasmodium falciparum in the host.An abundant literature has been dedicated to malaria's intra-host models, one can see for instance [1,5,12,14,15] and references therein.A review of some models has been given in [17].Most of these papers are, in the whole, related to the seminal Anderson, May and Gupta's model [1] which is the 3 dimensional system: where x is the concentration of uninfected erythrocytes in the blood, y is the concentration of infected erythrocytes, and m the concentration of free merozoites.Parameters µ x , µ y and µ m are the death rates of the uninfected erythrocytes, infected parasites and free merozoites respectively.The parameter β is the contact rate between erythrocytes and merozoites.Uninfected blood cells are recruited at a constant rate Λ from the bone marrow and have a natural life-expectancy of 1/µ x × days.Death of an infected erythrocyte results in the release of an average number of r merozoites.Free merozoites die or successfully invade a healthy erythrocyte.
For the epidemiology of malaria, we refer to [1,21].
As we can see in (1), the dynamics of the infected erythrocytes are described by one differential equation.However, experimental analysis has shown that the dynamics of Plasmodium falciparum might change according to its "age", i.e., the stage of its life cycle.Hence, a model, where the dynamics of the healthy cells are coupled with an age-structured model for the infected cells and the dynamics of free merozoites, has been considered [12,5]: One of the characteristics of Plasmodium falciparum is sequestration.To explain, let us give a brief description of the Plasmodium life cycle.The cycle begins when the parasite enters the human body through the bite of an infected mosquito, after which it migrates to and multiplies within the liver.Then the free forms (merozoites) resulting from this multiplication are able to invade the red blood cells (erythrocytes).Erythrocytes which are infected, mature during the erythrocytic cycle.At roughly the middle stage of trophozoite development (24 hours), molecules on the surface of infected erythrocytes link to receptors of endothelial cells.This bind has the effect of holding infected erythrocytes within vessels of organs (such as the brain), where they remain until the rupture of the erythrocyte and the release of merozoites.This period of attachment is called sequestration and during it, the infected erythrocytes are not detectable in the blood flow, they are "sequestered".Also it is widely accepted that antimalarial drugs act preferentially on different stage of parasite development [7,8].
In practice, to know a patient's stage of infection, the total parasite concentration n i=1 y i in the bloodstream is needed.However, only the peripheral infected erythrocytes (young parasites y 1 + y 2 + . . .y k , for some k < n), also called circulating, can be observed (seen on peripheral blood smears) and the other ones (sequestered: y k+1 , . . .y n ) are hidden in some organs like brain and heart, and cannot be observed.There is no clinical method of measuring the sequestered infected cells directly.
All parameters in (2) can be estimated by biological considerations, except the parameter β which is the Red Blood Cells (RBC)'s rate of infection.The rate of infection is generally an unknown in epidemiology models.
The estimation of sequestered parasite population has been a challenge for the biologist and modeler, with many authors [8,7,19] having studied this problem.In this paper, we propose a different and simple method to estimate the total parasite concentration from the measured circulating parasites and we show how the parameter β can be estimated.To this end, we use some tools from control theory.More precisely, we built an observer, i.e: an auxiliary dynamical system whose states tend exponentially to the states of the original system (2).Moreover, this auxiliary system does not depend on the badly known parameter β.The main idea here is to consider the infection term βxm as an unknown input and to use the theory of observers for systems with unknown inputs [10,11] to provide the dynamical estimates.
2. Dynamical estimation of the hidden parasitized erythrocytes.The exact number of stages of parasitized erythrocytes is unknown but according to [4] there are five main stages defined by simple morphology: young ring, old ring, trophozoite, early schizont and finally late schizont.So we suppose that the parasitized erythrocytes population within the host is divided in 5 different stages: y 1 , y 2 , y 3 , y 4 , y 5 .The two first stages correspond to the concentration of free circulating parasitized erythrocytes and the three last stages correspond to the sequestered ones.Let x be the concentration of healthy cells, and m the concentration of merozoites.The healthy cells x are produced by a constant recruitment Λ from the thymus and they become infected by an effective contact with a merozoite m.At the late stage of infected cells, the erythrocyte ruptures and releases r merozoites.The model that we consider has a seven dimensional state space, given by The different parameters of this model are defined as follows: • Λ : Recruitment of the healthy red blood cells (RBC).
• β : Rate of infection of RBC by merozoites.
• µ x : Natural death rate of healthy cells.
• µ i : Natural death rate of i-th stage of infected cells.
• γ i : Transition rate from i-th stage to (i + 1)-th stage of infected cells.
• r : Number of merozoites released by the late stage of infected cells.
• µ m : Natural death rate of merozoites.
As suggested in [4], we suppose that we can measure the circulating parasitaemia, i.e: y 1 + y 2 and we want to find an estimate of the sequestered parasitaemia, y 3 + y 4 + y 5 .Here, we use some classical notations from control theory.Let z(t) be the state at time t of the system (3), i.e., z = (x, y 1 , . . ., y 5 , m) ∈ R 7 .The measurable output of system (3) will be denoted by Y. Here, the measurable output corresponds to the concentration of the circulating parasitized erythrocytes, and hence we have Y(t) = y 1 (t) + y 2 (t).System (3) can then be written as: where z ∈ R 7 and Y ∈ R are the state vector and the measurable output, respectively.The matrices A, C, E and e 1 are defined as follows: , and C = (0, 1, 1, 0, 0, 0, 0).
The equation (4a) describes the dynamics of the system while equation (4b) gives information about what is measured.The problem we are interested in is: how to use the dynamical model (4a) together with the information provided by the measurements (4b) in order to get estimates of the state variables that cannot be measured?A solution to this problem can be provided by what is called in Control Theory an observer or estimator.An observer is an auxiliary dynamical system which uses the information Y(t) provided by the system and which produces an estimate ẑ(t) of the state z(t) such that the estimation error ẑ(t) − z(t) tends exponentially fast to zero as time t goes to infinity.An observer for system (4a-4b) will be of the form One has to find the functions f and g in such a way that the solutions ẑ(t) of ( 5) and the solutions z(t) of (4a) satisfy for all initial conditions (ẑ(0), z(0)): for some positive real number λ.
However, the construction of the observer (5) requires the knowledge of all the parameters of the model.This is not the case for our system since the parameter β is unknown.Therefore we are going to consider the term βxm appearing in the system (4a) as an unknown input and we shall use the notation: With this, system (4a) can be considered as a linear system with unknown input d(t) and can be written: Following the approach from [11], we shall construct an estimator (observer) for system (6).The idea for the construction will become evident after the proof of Theorem 2.1 below.
We define the matrix Ā by where I d denotes the identity matrix.For our system, we have It can be shown that the pair ( Ā, C) is detectable which means that there exists a column matrix such that the matrix Ā − LC is Hurwitz, i.e., all its eigenvalues have a negative real part.
We can now give our main result: we have all the ingredients to derive a method to dynamically estimate the state z(t) of our system.Theorem 2.1.An exponential observer for system (6) is given by This observer is inspired from [11].
Proof.We prove that the estimation error e(t) = ẑ(t)−z(t) converges exponentially towards zero.Indeed, we have This shows that e(t) → 0 with exponential convergence rate since all the eigenvalues of the matrix Ā − LC have negative real parts.
Remark 1.It is remarkable that the computation of the estimates ẑ(t) does not depend on the parameter β.Moreover, the estimation method proposed in Theorem 2.1 is valid for arbitrary number n of age-classes.
For our model with n = 5, we have , −µ m , and −(l 2 + l 3 ).Therefore the matrix Ā − LC is stable provided that l 2 + l 3 > 0. Since the eigenvalues of Ā − LC do not depend on l 1 nor on l i for i = 4, ..7, we simply choose l 1 = l 4 = l 5 = l 6 = l 7 = 0.The eigenvalues −(l 2 + l 3 ) and −(γ 1 + γ 2 + µ 2 ) come from the 3 × 3 upper left block, and it is somehow a stroke of luck that they are so simple algebraically.
Therefore one can get an estimate of the total parasite burden 5 i=1 y i (t).This dynamical estimate is given by ẑi (t).

3.
A method for the estimation of the rate of infection β.Let us consider the dynamic of y 1 , we have ẏ1 = βxm − (γ 1 + µ 1 ) y 1 .By the variation of constants method, we can write: x(s)m(s)e (γ1+µ1) (s−t) ds By replacing x(t), y 1 (t) and m(t) by their estimates x(t), ŷ1 and m(t) provided by the estimator (8) , we have: x(s) m(s)e (γ1+µ1) (s−t) ds Then ŷ1 (t)e (γ1+µ1) t − ŷ1 (t 0 )e (γ1+µ1) t0 = β t t0 x(s) m(s)e (γ1+µ1) s ds (9) Let t f be the time after which the measured peripheral parasitaemia Y becomes practically zero.We discretize the interval [t 0 , t f inal ] in [t i , t i+1 ] and we write the relation ( 9) for each interval [t i , t i+1 ] to obtain For each i, let ti , and x(s) m(s)e (γ1+µ1) s ds.
By numerical integration, we can compute, for each i, U i and V i .Thus, (10) can written as: U = βV where U and V are vectors of appropriate dimension with U i and V i as components of U and V respectively.Hence, by a linear regression, we can obtain β which is an estimate of β.

4.1.
Estimation of the numbers of sequestered infected erythrocytes.In this section we show how our method to estimate the total parasitaemia from the peripheral infected cells that are measured can be applied.By peripheral parasitaemia, we refer the density of parasites in the blood.The measured peripheral parasitaemia are taken from data collected when malaria was used as therapy for neurosyphilis by the US Public Health Service at the National Institutes of Health laboratories in Columbia, South Carolina and Milledgeville, Georgia [2].Patients were inoculated through mosquito bite or infected blood [2,3].Those data have been successfully used by many authors [2,3,6,18] for different purposes concerning the Plasmodium falciparum parasite.
To apply our total parasitaemia estimation method, we choose data of four patients named S1204, S1050, G221, and G54.These patients have been studied in [3] for the following reasons: they do not have any of the following features: • (i) treatment with drugs reported to affect gametocytes (chlorguanide, pyrimethamine, primaquine) or unknown drugs; • (ii) super-inoculation with any malaria parasite in the course of the infection; • (iii) parasitological status on the day before onset of curative treatment positive or uncertain; • (iv) gametocyte counts exceeding 100 µl (otherwise a reliable estimation of parameters would have been impossible with the current model); • (iv) inoculation of a rarely used strain.Those particular patients are relevant to our study because those patients have not taken treatment with drugs reported to affect gametocytes [3] and then satisfy the assumption of low mortality rates for peripheral erythrocytes that we make here.This assumption has also been used in [4].
It remains to pick the numerical values for different parameters.Among parameters in (3), some of them are known or at least widely accepted.However, the infection rate β is unknown or hardly known, so our method here allows us to study the system without explicitly using β.With many classes, the most difficult task is how to choose the transition rates [4,8] with the condition that the sum of the stage waiting times (life-span of parasitized erythrocytes) is 2 days.So, in this paper we use the transition and death rates from [4] for patients without fever, i.e: at 37 With these parameters, we will obtain 9101065 days, which is close to 2 days.In this paper, the unit of volume is micro-liter (µl) and the unit of time is day.
For the constant recruitment, we proceed as in [16], where the authors fixed the constant rate of erythrocytes' production at 5 × 10 6 120 cells µl −1 day −1 .This quantity is obtained from the parasite-free equilibrium state: 0 = Λ − µ x x 0 where x 0 is the fixed concentration of RBC in adult males which is around 5 × 10 6 RBCs per µl.

Merozoite's life-span ( 1 µ m
) is, according to [1], about 20 minutes, so µ m = 72 The number of merozoites released in the blood is between 8-32.So, we take here r = 16.
Recall that the choice l 2 and/or l 3 positive make the the matrix Ā − LC asymptotically stable.More precisely the spectrum of the matrix Ā − LC is: For the simulations we take l 2 = l 3 = 5.Since all the eigenvalues of the matrix Ā − LC are negative, the estimation error will tend exponentially to zero.For a given precision > 0, it is possible to compute a minimal time t satisfying | ŷi (t) − y i (t)| < ẑ(0) − z(0) for all t > t .In practice, to determine t , it is sufficient to run the estimator (8) with different initial conditions and then the estimates become reliable from the time when the various curves merge.An illustration for the total parasite burden is given in Figure 1.For the patient S1204, t ≈ 5 days.
Hereafter, we present the estimation results for the patient S1204 obtained by the estimator (8) using the malariatherapy data.The total parasitaemia estimates for the other patients are given in Figures 6,7,8.
The estimates of sequestered and total parasite population per mico-liter in patient S1204 are given in Table 1.4.2.Partial synchrony.One of the most important feature of malaria is the synchronous development of the parasite within the host, i.e. all of the infected RBCs (rings, trophozoites, shizontes) are about the same during the development of the cycle.This synchrony leads to simultaneous eruption of all mature schizontes, which explains the periodic fevers in human malaria [4,13].According to [20], synchronization in Plasmodium falciparum is only partial, i.e: parasiteamia often shows chaotic fluctuations on top of periodic behavior.This is confirmed by the malariatherapy patients considered in this paper.This can be seen on Fig. 6 for the patient G54: on the interval of (0,10), there are four main peaks of total parasiteamia.A high synchrony would suggest five peaks.The same remark is valid for the patients S1204 on Fig. 1 and S1050 on Fig. 7.For the patient G221, a closer look of the parasiteamia on the interval (0,10) on Fig. 9 shows also four peaks on that interval of time.

total parasitaemia per micro-liter t (days)
Another way to see the partial synchrony in the malariatherapy patients is to estimate the level of synchrony with the entropy formula [4] where p i is the proportion of stage i in the total infected population at time t: .
This formula suggests that if s(t) is close to zero, the population is highly synchronous (all the parasites are in the same age-class) and less synchronous otherwise (see, for instance, [4]).The following figures represent the level of synchrony in two of our four malariatherapy patients.These observation of low synchrony Synchrony Ratio would serve as an additional validation of the model since it is classically observed that Plasmodium falciparum malaria is less synchronous than other form of human malaria [13,20].Also, since synchrony and the mean age of parasites define how peripheral parasiteamia may reflect total parasiteamia [8], one can estimate the ratio of sequestered to circulating parasites, over time, as predicted by the model.When the patient is not cleared by parasiteamia, this ratio (Fig. 3) fluctuates over time.
For this particular patient, after the onset of the infection, the ration is greater than 1, i.e: there are more sequestered than circulating parasites.Another interesting feature that can be seen on Fig. 3 is that the ratio increases or decreases along with synchrony after the tenth day of infection.However at the onset of the infection the level of synchrony and ratio is more disparate (ratio less than 1 while the synchrony hits two main peaks).A possible explanation of this phenomenon is that at the beginning of the infection, the correlation between sequestered and circulating parasites depends more on the mean age of parasites than on synchrony [8].
4.3.Numerical estimation of β.Following Section 3 and by implementing the equation ( 10) on Maple, we obtain numerical estimated values of β for different patients.These estimates are given in Table 2.
Table 2.Estimated values for β based on the estimator (8) In the next section we provide an estimate using an aggregate model and compare with other values of β from the literature.
5. Aggregate model.The exact number of classes in the circulating or sequestered phases of parasites is not completely known.Hence, we divide the dynamics of parasite into two main classes, viz circulating parasites and sequestered parasites.Gravenor et al. [8] have considered this decomposition as well, but they consider only the dynamics of the infected cells.Figure 5 illustrates this decomposition.Therefore, we aggregate the different infected erythrocyte stages in two stages y 1 and y 2 describing the circulating parasites that are measured and the sequestered parasites respectively.We then obtain the following aggregated model: In our control theory framework, the system (11) could be written as: where: , n.At the half-way point of parasite t, the infected erythrocyte leaves the blood and binds to endothelium in asculature where the cycle is comghter parasites released at erythrore-enter the circulation and invade a rocyte.A measurement of P. falcisitaemia taken from a blood smear amples young parasites only.It is relate this measure to the total sity.In many cases a population of velops in synchrony.A low periphaemia may reflect low or high total mbers depending on the level of and the mean age of the parasite Peripheral parasitaemia therefore always to be a good correlate of meters.cult to form a reliable picture of the antimalarial therapy without knowing ur of the sequestered parasite populaantimalarial drugs are known to act y on different stages of parasite t, it is conceivable that a drug that ared parasites from the peripheral effect slower clearance of sequestered his is of particular importance since uestration is considered central to the f severe malaria.White et al. (1992) mathematical models can be used to terns of parasite sequestration.Grave-1998) presented a simple method for estimates of the level of sequestered m observed peripheral parasitaemia undergoing drug treatment.Here, we eneral approach to modelling the age P. falciparum that can be adapted to ticular data set that is to be analysed.a test of the approach using in vitro of parasites, apply the model to two nical data sets and use the model to investigate the relationship between erature and total parasite density.The same observer (7) is valid.The choice of parameters is slightly easier in this case.If we consider patients that are not taking drugs, the mortality of parasite is very low [4].Hence, we can consider, λ 1 = 1.03, λ 2 = 0.74 and µ 1 = 0.42, µ 2 = 0.08 .The unit remains days −1 .We choose these parameters on the basis of the precedent application.The same constraint 2 i=1 1 λi = 1.9101065 days is satisfied.The other parameters are standards and are the same as for the above application.With the choice L = (0, 5, 0, 0) T , we have successively:

Materials and Methods
Hence, for this model, the observer is given by: With this observer (which is independent of β), we can estimate all vector states and then the total parasite load y 1 + y 2 .Using this observer (13) with data corresponding to the four patients we obtain estimates of the total parasite load that are almost the same as those obtained with the first observer (8).
Following (10), we have implemented the observer using the data and then used numerical integration together with linear regression to estimate β.The obtained estimations of β for different patients with this method are given in Table 3.It can be noticed that the values (Table 3) obtained with the aggregated model are comparable to those (Table 2) obtained with the original system.Table 3.Estimated values for β based on the estimator (13).
It is worthwhile to compare these results with some others in the literature.Hetzel and Anderson [9] found, for rodents, β between 1.872 × 10 −6 (cells/µl) −1 day −1 and 7.329×10 −6 (cells/µl) −1 day −1 .Gravenor et al. in [6] have considered the model (1) and with a quasi steady-state assumption and by using clinical data, have estimated β = 3.4833 × 10 −6 (cells/µl) −1 day −1 .However as stated by the authors, those assumptions are simplifying.By exploiting clinical data of parasitaemia used in [6] in our method, we find β = 6.2734 × 10 −6 (cells/µl) −1 day −1 .According to [6], in order for the model to predict standard clinical parasitaemia, β would have to be less than 2 × 10 −6 (cells/µl) −1 day −1 (In fact β/µ m < 10 −3 in their units).Our results are consistent with their claim.In conclusion, one of our main objectives here is also to show how techniques from control theory are useful to estimate parameters of epidemic and intra-host models, and are therefore of wide applicability.

Figure 1 .Figure 2 .
Figure 1.Patient S1204: Estimates of the total burden with different initial conditions for the estimator (8): estimates become accurate for t > 5 days.

Figure 3 . 4 Figure 4 .
Figure 3. Patient G54: Degree of population synchrony (solid line) and ratio of sequestered to circulating parasites (dashed line).In this case the maximum synchrony is -1.129 while the minimum is -1.571.

Fig. 1 .
Fig.1.The mathematical model of the life cycle is based on a finite number of compartments (here depicted as circles), each representing an equal duration of development time.In the above example, there are eight Compartments and since the parasite life cycle is 48 hr, each compartment represents 6 hr.Parasites can often be aged on appearance.Four commonly used morphological stages are young rings, late rings/young trophozoites, old trophozoites and schizonts/segmenters.In this model, these are each represented by two compartments (though more can be used).In an infected individual, parasites in approximately the first half of the cycle (young rings-young trophozoites) circulate freely and can be seen in the peripheral blood, while all other parasites sequester in the deep vasculature (old trophozoites-segmenters) and cannot be detected.

Figure 5 .
Figure 5. Decomposition of Plasmodium falciparum cycle into two classes.