Quantification of Healthy Red Blood Cell Removal and Preferential Invasion of Reticulocytes in Macaca mulatta during Plasmodium cynomolgi Infection

TL;DR

This study develops and calibrates a simplified mathematical model to quantify red blood cell removal and reticulocyte invasion during malaria infection in macaques, revealing host erythropoietic response prior to blood stage.

Contribution

The paper introduces a simplified ODE model derived from a PDE framework and calibrates it with experimental data to quantify key malaria dynamics in macaques.

Findings

Quantified healthy red blood cell removal during infection.

Identified preferential invasion of reticulocytes.

Revealed host erythropoietic response before blood stage.

Abstract

We derived an ordinary differential equation model to capture the disease dynamics during blood-stage malaria. The model was directly derived from an earlier age-structured partial differential equation model. The original model was simplified due to experimental constraints. Here we calibrated the simplified model with experimental data using a multiple objective genetic algorithm. Through the calibration process, we quantified the removal of healthy red blood cells and the the preferential infection of reticulocytes during \textit{Plamodium cynomolgi} infection of \textit{Macaca mulatta}. The calibration of our model also revealed the existence of host erythropoietic response prior to blood stage infection.

Figures

Tables

Table 1. Table 1: Enrichment of GO Gene Sets

Positive Enriched GO Gene Sets	Negative Enriched GO Gene Sets
INNATE IMMUNE RESPONSE	NCRNA PROCESSING
IMMUNE EFFECTOR PROCESS	NCRNA METABOLIC PROCESS
ACTIN FILAMENT BASED PROCESS	RIBOSOME BIOGENESIS
REGULATION OF BODY FLUID LEVELS	RRNA METABOLIC PROCESS
WOUND HEALING	RIBONUCLEOPROTEIN COMPLEX BIOGENESIS
RESPONSE TO BACTERIUM	TRANSLATIONAL INITIATION
DEFENSE RESPONSE TO OTHER ORGANISM	AMIDE BIOSYNTHETIC PROCESS
INFLAMMATORY RESPONSE	PEPTIDE METABOLIC PROCESS
ACTIVATION OF IMMUNE RESPONSE	MULTI ORGANISM METABOLIC PROCESS
REGULATION OF RESPONSE TO WOUNDING	NONSENSE MEDIATED DECAY
HEMOSTASIS	PROTEIN LOCALIZATION
REGULATED EXOCYTOSIS	PROTEIN LOCALIZATION TO ENDOPLASMIC RETICULUM
SECRETION	PROTEIN TARGETING TO MEMBRANE
RESPONSE TO WOUNDING	TRNA METABOLIC PROCESS
CELLULAR RESPONSE TO NITROGEN COMPOUND	RNA CATABOLIC PROCESS
ENDOCYTOSIS	ORGANIC CYCLIC COMPOUND CATABOLIC PROCESS
RESPONSE TO VIRUS	TRNA PROCESSING
CELL CELL ADHESION	MITOCHONDRIAL TRANSLATION
REGULATION OF INFLAMMATORY RESPONSE	TRANSLATIONAL ELONGATION

Table 2. Table 2: Enrichment of Pathways

Positive Enriched Pathways	Negative Enriched Pathways
SYSTEMIC LUPUS ERYTHEMATOSUS	3 UTR MEDIATED TRANSLATIONAL REGULATION
HEMOSTASIS	TRANSLATION
CYTOKINE SIGNALING IN IMMUNE SYSTEM	PEPTIDE CHAIN ELONGATION
CHEMOKINE SIGNALING PATHWAY	RIBOSOME
INTERFERON ALPHA BETA SIGNALING	INFLUENZA VIRAL REPLICATION
PID PDGFRB PATHWAY	PROTEIN TARGETING TO MEMBRANE
REGULATION OF ACTIN CYTOSKELETON	NONSENSE MEDIATED DECAY
PLATELET ACTIVATION SIGNALING AND AGGREGATION	INFLUENZA LIFE CYCLE
AMYLOIDS	METABOLISM OF PROTEINS
PID VEGFR1 2 PATHWAY	METABOLISM OF RNA
RNA POL I PROMOTER OPENING	FORMATION OF THE TERNARY COMPLEX
NABA MATRISOME ASSOCIATED	ACTIVATION OF THE MRNA BINDING TO 43S
RESPONSE TO ELEVATED PLATELET CYTOSOLIC CA2	METABOLISM OF MRNA
NABA SECRETED FACTORS	MITOCHONDRIAL PROTEIN IMPORT
INTERFERON SIGNALING	RNA POLYMERASE
RHO PATHWAY	RNA POL III TRANSCRIPTION INITIATION
LEUKOCYTE TRANSENDOTHELIAL MIGRATION	RESPIRATORY ELECTRON TRANSPORT ATP SYNTHESIS
FOCAL ADHESION	PYRIMIDINE METABOLISM
JAK STAT SIGNALING PATHWAY	TCA CYCLE AND RESPIRATORY ELECTRON TRANSPORT

Equations

$$\begin{array}{rl}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial a}=& -\left(\sum _{i=1}^T{w}_i(t){\theta }_i+\sum _{i=1}^Q{s}_i(t){\psi }_i+h(a)\right)u(a,t)\\ & -\gamma \kappa v({\alpha }_{max},t)p(u(a,t)),\\ \frac{\partial v}{\partial t}+V\frac{\partial v}{\partial \alpha }=& -\left(\sum _{i=1}^T{w}_i(t){\varphi }_i+\sum _{i=1}^Q{s}_i(t)b_i(t)\right)v(\alpha ,t),\\ \frac{d{w}_i}{dt}=& o_i(t)-{\beta }_i{w}_i(t),\\ \frac{d{s}_i}{dt}=& l_i(t)-{\delta }_i{s}_i(t),\end{array}$$

$$\begin{array}{rl}u(a,0)=& g(a),\\ u(0,t)=& f(t,\phi (t)),\\ v(\alpha ,0)=& c(\alpha ),\\ v(0,t)=& \gamma \kappa {\int }_{a_0}^{a_{max}}v({\alpha }_{max},t)p(u(a,t))da.\end{array}$$

$$\begin{array}{rl}\frac{\partial u}{\partial t}+\frac{\partial u}{\partial a}=& -h(a)u(a,t)-\gamma \kappa v({\alpha }_{max},t)p(u(a,t)),\\ \frac{\partial v}{\partial t}+V\frac{\partial v}{\partial \alpha }=& 0.\end{array}$$

$$V(t)={\int }_0^{{\alpha }_{max}}v(\alpha ,t)d\alpha ,$$

$$V(0)={\alpha }_{max}C.$$

$$\begin{array}{rl}{\int }_0^{{\alpha }_{max}}v(\alpha ,t_2)d\alpha =& {\int }_0^{{\alpha }_{max}}v(\alpha ,t_1)d\alpha ,\\ & +{\int }_{t_1}^{t_2}\gamma \kappa {\int }_{a_0}^{a_{max}}v({\alpha }_{max},t)p(u(a,t))dadt,\\ & -{\int }_{t_1}^{t_2}v({\alpha }_{max},t)dt.\end{array}$$

$$\begin{array}{rl}{\int }_0^{{\alpha }_{max}}v(\alpha ,t_2)d\alpha =& {\int }_0^{{\alpha }_{max}}v(\alpha ,t_1)d\alpha ,\\ & +{\int }_{t_1}^{t_2}\gamma \kappa v({\alpha }_{max},t)dt,\\ & -{\int }_{t_1}^{t_2}v({\alpha }_{max},t)dt.\end{array}$$

$$V(t_2)-V(t_1)=(\gamma \kappa -1){\int }_{t_1}^{t_2}v({\alpha }_{max},t)dt.$$

$${\int }_{t_1}^{t_2}\dot{V}dt=(\gamma \kappa -1){\int }_{t_1}^{t_2}v({\alpha }_{max},t)dt$$

$$\dot{V}=(\gamma \kappa -1){\int }_{a_0}^{a_{max}}v({\alpha }_{max},t)p(u(a,t))da.$$

$$R(t)={\int }_0^{a_1}u(a,t)da,$$

$$U(t)={\int }_{a_1}^{a_{max}}u(a,t)da,$$

$$R(t)+U(t)={\int }_0^{a_{max}}u(a,t)da.$$

$$\dot{V}=(\gamma \kappa -1)v({\alpha }_{max},t)\left(\frac{R}{R+U}+\frac{U}{R+U}\right).$$

$$v({\alpha }_{max},t)=C_{2}V.$$

$$\dot{V}=(\gamma \kappa -1)C_{2}V\left(\frac{R}{R+U}+\frac{U}{R+U}\right).$$

$$\dot{V}=(\gamma {\kappa }_1-1)C_{2}V\frac{R}{R+U}+(\gamma {\kappa }_2-1)C_{2}V\frac{U}{R+U}.$$

$$x_{3}N=(\gamma {\kappa }_1-1)C_2,$$

$$x_{6}N=(\gamma {\kappa }_2-1)C_2.$$

$$\dot{V}=x_{3}N\frac{UV}{R+U}+x_{6}N\frac{RV}{R+U}.$$

$$\begin{array}{rl}{\int }_0^{a_1}u(a,t_2)da=& {\int }_0^{a_1}u(a,t_1)da\\ & +{\int }_{t_1}^{t_2}{x}_4{e}^{x_5(T_0-T)}dt\\ & -{\int }_{t_1}^{t_2}u(a_1,t)dt\\ & -{\int }_{t_1}^{t_2}{x}_{6}N\frac{RV}{R+U}.\end{array}$$

$$R(t_2)-R(t_1)={\int }_{t_1}^{t_2}{x}_4{e}^{x_5(T_0-T(t))}-u(a_1,t)-x_{6}N\frac{RV}{R+U}dt.$$

$$\dot{R}=x_4{e}^{x_5(T_0-T)}-u(a_1,t)-x_{6}N\frac{RV}{R+U}.$$

$$\dot{R}=x_4{e}^{x_5(T_0-T)}-x_{1}R-x_{6}N\frac{RI}{R+U}.$$

$$\begin{array}{rl}{\int }_{a_1}^{a_{max}}u(a,t_2)da=& {\int }_{a_1}^{a_{max}}u(a,t_1)da\\ & +{\int }_{t_1}^{t_2}u(a_1,t)dt\\ & -{\int }_{t_1}^{t_2}u(a_{max},t)dt.\\ & -{\int }_{t_1}^{t_2}{\int }_{a_1}^{a_{max}}h(a)u(a,t)dadt\\ & -{\int }_{t_1}^{t_2}{x}_{3}N\frac{UV}{R+U}.\end{array}$$

$$U(t_2)-U(t_1)={\int }_{t_1}^{t_2}u(a_1,t)-x_{2}U-x_{6}N\frac{RI}{U+R}dt.$$

$$\dot{U}=x_{1}R-x_{2}U-x_{6}N\frac{RI}{U+R}.$$

$$\begin{array}{rl}\dot{U}& =x_{1}R-x_{2}U-x_{3}N\frac{UI}{T},\\ \dot{R}& =x_4{e}^{x_5(T_0-T)}-x_{1}R-x_{6}N\frac{RI}{T},\\ \dot{I}& =x_{3}N\frac{UI}{U+R}+x_{6}N\frac{RI}{U+R}.\end{array}$$

$$x_4=x_1{R}_0=x_2{U}_0.$$

$$AP{E}_{retic}=\frac{100}{n}\sum _{i=1}^n\frac{r_i-R_i}{R_i},$$

Full text

1 Introduction

The goal of this article is to quantify the removal of healthy red blood cells and preferential invasion of reticulocytes during malaria infection, and find their associations with molecular phenomena through a modeling approach. The disease malaria is caused by Plasmodium parasites. Out of the five human Plasmodium species capable of causing malaria, Plasmodium falciparum and Plasmodium vivax account for the majority of human malaria infections. P. falciparum is responsible for the majority of malaria-related mortality and is most prevalent in sub-Saharan Africa [1].

The Plasmodium life cycle is comprised of several stages. The infection process in humans starts with the injection of sporozoites by mosquitoes into the skin of the host. This is followed by the liver stage, during which the inoculated sporozoites grow and multiply asexually within hepatocytes for 1-2 weeks to produce merozoites. The newly produced merozoites emerge from the liver and enter the blood stream. The blood-stage infection starts immediately after the hepatic stage. The parasite’s blood-stage infection in both human and non-human primates generally has a regular cycle of 24 to 72 hours depending on the species of the Plasmodium parasite [2, 3]. The parasites invade healthy red blood cells (RBCs) and replicate asexually, remodeling and ultimately destroying the RBCs in the process. The destruction of RBCs during blood-stage malaria infection sometimes results in severe anemia, which is one of the major complications of malaria and a leading cause of mortality.

During blood-stage malaria infection, the majority of red blood cell loss has been attributed to the removal of healthy red blood cells [4, 5, 6]. The hemodynamic model described by Yan et al [7] characterizes this phenomenon through a mechanistic model taking into account the interaction between healthy RBCs, infected RBCs, cells form the innate immune system, and cells from the adaptive immune system. This model captures the clinical outcomes of a malaria infection: resistance to the disease, susceptibility, and resilience (chronic infection with mild symptoms). That model has been simplified to account for experimental constraints, and it has been calibrated in this article with the data of Macaca mulatta infected with P. cynomolgi described by Joyner et al [8]. Plasmodium cynomolgi is a malaria parasite that infects old world monkeys; it is physiologically and evolutionarily similar to P. vivax [9, 10].

This paper is organized as follows: Section 2 provides the description of the experimental data. Section 3 contains the detailed derivation and description of our hemodynamic model. Section 4 describes the model calibration process. Section 5 provides an overview of the calibration results. Section 6 integrates our modeling result with transcriptomic data collected during the experiment. Section 7 describes an adjusted model based on our calibration results. Section 8 offers some conclusions, and discusses the biological significance of our results.

2 Experimental Description

The experimental design is shown in Fig 1. The infected red blood cell, total red blood cell, and reticulocyte concentration were measured daily over the entire experiment. The four subjects, S1, S2, M1 and M2 are referred as Subject 1 $\sim$ 4 within this paper. Each red diamond corresponds to a specific time point where transcriptome, cytokine and immune profiling was conducted. Time point 2 (TP2) corresponds to the acute primary infection when parasite counts peaked. A detailed description of the experiment that collected this data can be found in Joyner et al. [8].

3 Matching a Theoretical Model with Experimental Constraints

The mechanistic partial differential equation model described by Yan et al. [7], has the form presented below; the detailed definition of variables, parameters and functions is presented in that paper, and also for completeness in the Appendix.

[TABLE]

subject to the following initial and boundary conditions:

[TABLE]

With the experimental constraints, it was not possible to measure the delay in the erythropoietic response, and hence it was removed from this model. Thus the boundary condition $u(0,t)=x_{4}e^{x_{5}(T_{0}-T(t))}$, where $x_{4}$ is the baseline production of red blood cells and $x_{5}$ is a parameter controlling the speed of erythropoiesis response. $T_{0}$ is the sum of the steady-state concentration of reticulocytes and mature red blood cell and $T$ is the sum of the concentration of reticulocytes and mature red blood cell at time $t$.

Under the explicit assumption, that the immune cell populations do not interact with the infected red blood cell population ($v$) and healthy red blood cell population ($u$) ($\phi_{i}=0$, $\theta_{i}=0$, $\psi_{i}=0$, for all $i$ and $b_{i}(t)=0$ for all $t$), the system reduces to:

[TABLE]

Furthermore, under the explicit assumption that the initial age distribution of v is uniform, $v(\alpha,0)=C$. Let

[TABLE]

where $V(t)$ denote the total amount of iRBC at time $t$. Then,

[TABLE]

Additionally, let $t_{2}=t_{1}+\epsilon,\epsilon>0$, then:

[TABLE]

Because $\int_{a_{0}}^{a_{max}}p(u(a,t))=1$ for all $t$, the system further simplifies to:

[TABLE]

Substituting $V$ into (2), and rearranging the equation:

[TABLE]

Thus,

[TABLE]

and

[TABLE]

Because the age distribution of red blood cells are unknown throughout the experiment, and only reticulocyte and mature red blood cell populations are measured, we define $a_{1}$ as the precise time point where reticulocytes become mature red blood cells. We further define the total population of reticulocytes at time point $t$ as $R(t)$ and the total population of mature red blood cells at time point $t$ as $U(t)$. $R(t)$ and $U(t)$ have the following form:

[TABLE]

and

[TABLE]

thus

[TABLE]

Therefore, (3) can be rewritten as:

[TABLE]

Under the assumption that on average a small fixed percentage of infected red blood cells are bursting, then

[TABLE]

Which means:

[TABLE]

Additionally, knowing that the infection rate of reticulocyte is different from mature red blood cells, we can modify the system to reflect this fact:

[TABLE]

Let

[TABLE]

and

[TABLE]

Then

[TABLE]

The derivative of $R(t)$ can be derived similarly. Let:

[TABLE]

Because the much shorter life span of reticulocytes in comparison to mature red blood cells [11] and the fact that the survival rate of reticulocytes to be close to 1 [12], the hazard function, h(a), of reticulocyte is ignored. Substituting R, and rearrange the equations, we obtain:

[TABLE]

Thus,

[TABLE]

Under the assumption that on average a small fixed percentage of reticulocyte are aging into mature red blood cell at any given time,

[TABLE]

The derivative of $U(t)$ is also derived similarly. let:

[TABLE]

Assuming that on average, a fixed percentage of red blood cells are removed due to random chance, substituting U then the equation becomes:

[TABLE]

Thus:

[TABLE]

In conclusion, under the following explicit assumptions:

• •

The change of infected red blood cell population and healthy red blood cell population is independent of the immune cell populations,

• •

the erythropoiesis response does not have a delay,

• •

on average, a fixed percentage of iRBCs are bursting at any given moment and a fixed percentage of reticulocytes are transitioning into mature red blood cells,

the original PDE model (1) is simplified to the following system of ordinary differential equations that captures the change in the total infected red blood cell, red blood cell, and reticulocyte population:

[TABLE]

Where U denote the concentration of healthy mature red blood cells (RBCs), R denotes healthy reticulocytes (RT) and I denote the concentration of infected red blood cells (iRBCs). Let $T=U+R$ and $T_{0}=R_{0}+U_{0}$ where $R_{0}$ and $U_{0}$ denote the steady state concentration of RTs and RBCs in the absence of malaria infection. N denote the average merozoites produced per infected red blood cell. $x_{1}R$ describes the aging of RTs to become RBCs. $x_{2}U$ denote the random removal of mature red blood cells. $x_{3}$ and $x_{6}$ denote the infection success rate of RBCs and RTs respectively. $x_{4}$ denotes the baseline production of RTs and the term $e^{x_{5}(T_{0}-T)}$ describes the host erythropoiesis response.

4 Parameter Estimation

Under the assumption of steady state, where $R_{0}$ is the steady state concentration of RTs and $U_{0}$ is the steady state concentration of RBCs, the following equality is established:

[TABLE]

The average number of merozoite produced is set to be 20. Fonseca et al. [12] estimated that the baseline production of RTs in M mulatta, $x_{4}$, is $192,00$ cells per day [12]. $R_{0}$ and $U_{0}$ are estimated for each subjects using the average of the red blood cell and reticulocytes counts during the first ten days of the experiment where no iRBCs were detected. The lower bound and upper bound of the other three parameters, $x_{3}$, $x_{6}$ and $x_{5}$ were also estimated. All three parameters are positive. The upper bound of $x_{5}$ was set to allow a maximum of 2-fold increase in RTs production in each subject. Upper bound of $x_{3}$ and $x_{6}$ was set to be $0.05$ and $1$ respectively. When $x_{3}>0.05$ and/or $x_{6}>1$, the system becomes numerically unstable to solve due to stiffness using odesolver45 in MATLAB environment. The upper bound of $x_{3}$ and $x_{6}$ also ensures that the $\dot{I}<0.5NI$ throughout the simulation, meaning that the number of infected red blood cells produced over a two day period can not exceed the number of merozoites produced.

A multiple objective genetic algorithm [13] was used to estimate the parameters ($x_{3},x_{6},x_{5}$) for each subjects using the estimated lower and upper bound of each parameter. The model was fitted to the cellular data with starting date corresponding to the first appearance of iRBCs and ending date corresponding to time point 2. The two objectives minimized were the average percent error (APE) of predicted RTs and iRBCs concentration. APE have the following form for a specific variable such as RTs:

[TABLE]

where $r_{i}$ is the observed RTs concentration at time point $i$ and $R_{i}$ is the model predicted RTs concentration. Contrary to single objective minimization, a Pareto front for both objectives is estimated during each application of the genetic algorithm. The Pareto front refers to a set of possible values of both object function such that the decrease in one objective function necessitates the increase of the other [13]. Each application of a genetic algorithm is terminated when the change in the estimated Pareto front is less than the predefined tolerance ($0.0001$). Due to the stochastic nature of the genetic algorithm, 1,000 runs were applied to estimate the parameters for each subject. $>99\%$ of the application of genetic algorithm terminated due to the convergence of Pareto front, the rest did not converge during the maximum allowed run time of $60$ seconds.

5 Results

5.1 Result Overview

The Pareto front of iRBC APE and RT APE for each subject are shown in Figure 2. For subject 1, 2 and 4, there exist a sub-region on the Pareto front such that the iRBC APE and RT APE have a negative linear relationship, which means that $iRTAPE+RTAPE$ is close to constant. The Pareto front estimated for subject 3 does not contain such a region, indicating a lack of model fit for that specific subject. The iRBC APE for all four subjects has a range of $(45\%\sim 90\%)$ where the RT APE have a range of ($15\%\sim 50\%$). The simulation of the top 100 model with the lowest iRBCs APE are shown in Figure (3, 4, 5, 6).

5.2 Quantification of The Preferential Infection of Reticulocyte by P. cynomolgi

To quantify the preferential infection of RTs by P cynomolgi, the ratio of $x_{6}$ to $x_{3}$ in the top 100 model with the lowest iRBCs APE for each subject was calculated. The ratio $\frac{x_{6}}{x_{3}}$ is interpreted as the likelihood of P cynomolgi infecting RTs over infecting RBCs. The distribution of $\frac{x_{6}}{x_{3}}$ is shown in Figure 7. The means of the four distribution are 13, 17, 28 and 18 respectively.

5.3 Quantification of Removal of Healthy Red Blood Cell (hRBC)

The loss of healthy RBC (hRBC) is also calculated based on the 100 model prediction with the lowest iRBC APE. Total Loss of hRBC and Ratio of loss of hRBC and parasitized RBC are shown in Figure 8. On average, 1.5 million hRBCs per $\mu$l are cleared by the host through out the onset of the disease till acute primary infection. Furthermore, our model predicted that for each RBC parasitized, the host removes 3 hRBC. The quantification of the removal of hRBCs over the entire infection is shown in Figure 9. The model fit of RBC and retic Population after adjusting for removal of hRBCs are shown in Figure 10.

6 Correlation and Enrichment Analysis of Rate of hRBC Removal

The rate of hRBC removal was estimated for all four subjects (Fig 11) based on the estimated loss of hRBCs (Fig 9). Pearson’s correlation between hRBC removal rate and their corresponding transcript, immune cell and cytokine abundance was calculated. The p-value of the correlation was adjusted for false discovery rate and the transcript, immune cell and cytokine exhibiting significant correlation with hRBC removal rate (q-val $\leq 0.05$) are shown in Fig 12. Several pro-inflammatory ctokines such as IL6 and IL1B are positively correlated with the rate of hRBC removal. On the other hand, CD8 T cell population displayed significant negative correlation with hRBC removal rate.

Only seven genes: MYO3B, GAN, DNAJB4, TRIM45, TMEM150A, IL23R and BMF have shown significant correlation with the rate of hRBC removal. To fully explore the association of transcriptome change and rate of hRBC removal, Gene Set Enrichment Analysis [14] was applied to the correlation ranked gene lists. The most significantly enriched GO gene sets and pathways are shown in table 1 and table 2. Innate immune related gene sets and pathways are positively correlated with hRBC removal rate, which corresponds to the positive correlation between hRBC removal and inflammatory cytokines. On the other hand, transcripts that have negative correlation with hRBC removal are enriched in RNA and protein processing related pathways.

7 Empirical Model Adjustment

The simulation of the top 100 model with the lowest iRBC APE for subject 1,2 and 3 (Figure 3, 4, 5, 6) failed to capture the increase of RBC population during the early stages of malaria infection. Considering that the difference between RBC population and steady state RBC population is minimal during the early stage of malaria infection, this discrepancy suggest that RBCs are released during the early stage of blood stage malaria. Additionally, our correlation analysis identified several cell types, cytokines and transcripts (q-value $\leq 0.05$) that are linearly correlated with hRBC removal rate. This finding suggest the possibility of using these entities to predict hRBC removal rate. Taking these findings into consideration, our original model can be further expanded to:

[TABLE]

is linearly explained by a combination of the abundances of inflammatory cytokines. Where each $E_{i}$ denote the abundance of a cytokine exhibiting linear correlation with the rate of hRBC removal (Fig: 12) and each $\gamma_{i}$ quantifies the linear dependency. The term $g(M)$ denote the increase of RBC population during the early stage of malaria infection, that is dependent on some molecular quantity $M$, possibly pathogen related. The estimated form of $g(M)$ are shown in Figure 13 and the post-hoc fitted model is shown in Figure 14.

8 Discussion

In this chapter, a simplified Ordinary Differential Equation model was fitted to the time series data collected from the MaHPIC experiment. Our simplified model was derived from the earlier PDE system under three explicit assumptions. The simplification of the original model reduced the unknown parameters to 3 parameters with well-defined lower and upper bounds.

The model was fitted using multiple objective genetic algorithm and the Pareto front estimated (Fig: 2) demonstrate that the simplified model describes the dynamics of RTs much better than that of iRBCs as shown by the 2-fold difference in the range of RT APE and iRBC APE. The iRBC APE for all four subjects has a range of $(45\%\sim 90\%)$ where as the RT APE have a range of ($15\%\sim 20\%$). This result suggests that the immune function omitted in our model have a large impact on parasite population during the primary infection.

The best-fitted models were used to quantify the preferential infection of RTs by P. cynomolgi. Our model predicts that, on average, P. cynomolgi merozoites are 20 times more likely to infect RTs than RBC (Fig 7). The experimental verification of our prediction is difficult due to the lack of an in vitro system to study P. cynomolgi, but our prediction provides a general range of preferential infection of RTs, which can be utilized in future experimental designs and modeling studies.

Additionally, in all the best fit models (Fig 3, 4, 5, 6), our prediction of RBC concentration and RT concentration are higher than observed. This over estimation of RBC and RT concentration constitutes the majority of APE. Furthermore, our iRBC concentration predictions are on average lower than the observed value. We utilized this fact to estimate the lower bound of the amount of hRBCs removed by the host. Our model estimates that at least $50\%\sim 80\%$ of RBC loss during the primary infection of P. cynomolgi are due to the removal of hRBCs 8. Finally, our model shows that the speed at which hRBCs are removed increases throughout the infection.

Using the estimated hRBC removal rate, a correlation study was conducted to identify transcript, immune cell and cytokine abundance that has significant correlation with hRBC removal rate. Interestingly, innate immune related gene sets (Interferon Response, Rho Pathway and JAK-STAT pathway et al.) along with pro-inflammatory cytokines (IL-1B, IL-6 et al.) displayed significant positive correlation with hRBC removal. The association of severe malaria anemia with pro-inflammatory response has long been studied [15]; our analysis provides a list of possible cytokine biomarkers for the estimation of host clearance of hRBCs. Due to the experimental constraint where cytokine profiling was only conducted during peak parasitemia, our analysis lacks the resolution to provide a mechanistic explanation for the correlation between pro-inflammatory cytokines and host removal of hRBC. Furthermore, the observation that inflammation related genes and cytokines are differentially up-regulated in the severe subjects (Subject 1 and 2), along with the fact that severe subjects have a higher rate of hRBC removal during peak parastemia, suggest the possible role of inflammation associated hRBC clearance and clinical severity.

In conclusion, we have demonstrated that a simplified model with only three unknown parameters can be used to predict RT concentration with an APE of ($15\%\sim 20\%$). Despite the model’s relatively poor performance at predicting iRBC dynamics (APE $50\%\sim 70\%$), it can be used to estimate the preferential infection of RTs and hRBC removal during malaria infection. The estimation of the hRBC removal rate using our model along with the down-stream enrichment analysis reveals associations of hRBC removal and both the inflammatory response and CD 8 T-cell response. Application of this model to more time series data sets of malaria infection involving a variety of malaria species is necessary to validate our findings.

APPENDIX

The parameters, functions, and variables described by Yan et al. [7] are presented here for completeness. The reader is encouraged to read that article for a detailed description of each element.

[TABLE]

[TABLE]