Computational synchronization of microarray data with application to Plasmodium falciparum

Abstract

Background

Microarrays are widely used to investigate the blood stage of Plasmodium falciparum infection. Starting with synchronized cells, gene expression levels are continually measured over the 48-hour intra-erythrocytic cycle (IDC). However, the cell population gradually loses synchrony during the experiment. As a result, the microarray measurements are blurred. In this paper, we propose a generalized deconvolution approach to reconstruct the intrinsic expression pattern, and apply it to P. falciparum IDC microarray data.

Methods

We develop a statistical model for the decay of synchrony among cells, and reconstruct the expression pattern through statistical inference. The proposed method can handle microarray measurements with noise and missing data. The original gene expression patterns become more apparent in the reconstructed profiles, making it easier to analyze and interpret the data. We hypothesize that reconstructed gene expression patterns represent better temporally resolved expression profiles that can be probabilistically modeled to match changes in expression level to IDC transitions. In particular, we identify transcriptionally regulated protein kinases putatively involved in regulating the P. falciparum IDC.

Results

By analyzing publicly available microarray data sets for the P. falciparum IDC, protein kinases are ranked in terms of their likelihood to be involved in regulating transitions between the ring, trophozoite and schizont developmental stages of the P. falciparum IDC. In our theoretical framework, a few protein kinases have high probability rankings, and could potentially be involved in regulating these developmental transitions.

Conclusions

This study proposes a new methodology for extracting intrinsic expression patterns from microarray data. By applying this method to P. falciparum microarray data, several protein kinases are predicted to play a significant role in the P. falciparum IDC. Earlier experiments have indeed confirmed that several of these kinases are involved in this process. Overall, these results indicate that further functional analysis of these additional putative protein kinases may reveal new insights into how the P. falciparum IDC is regulated.

Introduction

Approximately 40% of the global population is at risk for contracting malaria, and an estimated 780,000 people die annually from this disease [1]. Human malaria is caused by five Plasmodium species, of which P. falciparum is responsible for the majority of human fatalities. The disease is transmitted when an infected mosquito bites a person, and injects sporozoites that migrate to and develop in the liver before merozoites are released into the bloodstream and invade red blood cells (RBCs) [2]. Within the RBC, P. falciparum undergoes a well-defined developmental cycle (IDC) during a 48-hour period that is characterized by three main stages, namely: rings, trophozoites and schizonts [2]. Schizont-infected RBCs rupture at the end of this cycle to release merozoites that can invade RBCs and reestablish a new IDC. In addition to the morphologic changes that characterize parasite development during the IDC, changes in gene expression accompany [3, 4] and most likely drive this developmental program.

Gene expression during the 48-hour IDC has been densely profiled at 1-hour intervals using microarray technology in an effort to understand how overall gene expression patterns help shape blood stage parasite biology [3, 4]. These studies revealed that the levels of many transcripts are reproducibly high or low at characteristic times within the IDC. Genes involved in key biological processes most relevant to a given IDC stage are generally coordinately up- and down- regulated. In fact, a 'just-in-time' model to describe transcriptionally regulated gene expression in P. falciparum has been proposed [4, 5]. Here, a transcript's level is proposed to peak just prior to when its encoded protein product is most critically needed. Regulation of general biological processes such as metabolism, DNA synthesis, protein turnover and red blood cell invasion, for example, is well-described by this model [3, 4]. While the full complement of parasite proteins controlling IDC progression has not been identified, the just-in-time principle could be useful for identifying key, transcriptionally regulated proteins that play important roles in regulating this process.

Protein kinases represent one such protein class, and have previously been implicated in regulating various aspects of the cell cycle and development in P. falciparum [6–8]. The genome encodes a predicted 85 [9] or 99 [10] protein kinases, which include an expanded and divergent FIKK family, of which there are 20 members [9, 11]. In Plasmodium spp., several protein kinases have been shown to be important in regulating blood stage biology, including parasite egress from infected RBCs [6–8]. Recently, a large-scale knockout screening effort identified 36 out of the 65 protein kinases evaluated (no FIKKs included) as likely essential to the development of P. falciparum blood stage parasites [12]. Overall, these studies highlight the important contribution of the protein kinases to regulating critical aspects of P. falciparum biology during the IDC.

Therefore, we have been interested in addressing whether computationally applying the just-in-time principle to publicly available microarray data for P. falciparum is a reasonably efficient strategy for identifying protein kinases that are important to transition through the IDC. Such an approach, if successful, could prioritize protein candidates for more exhaustive experimental analysis. Two fundamental assumptions underlie our analytical framework, namely: (1) an important subset of protein kinases regulating transition through the IDC is transcriptionally regulated, and this is captured in the publicly available microarray data; and (2) increases in protein kinase transcript levels predictably precede peak protein synthesis, the latter defining the time at which a given protein kinase plays its critical role in IDC progression. Our approach requires explicitly addressing the confounding factor of decaying synchrony of parasite cultures in order to improve the reliability of predictions. For microarray experiments, parasite cultures are initially synchronized. However, these cultures gradually lose synchrony over the experimental course [3]. Consequently, gene expression levels at discrete time points for individual parasites are not directly inferred from the observed microarray data, which reflect an ensemble of increasingly asynchronous parasite transcription profiles. To address this issue, we introduce a deconvolution approach for reconstructing the intrinsic gene expression pattern of individual parasites from microarray data. This work was partly presented in our previous paper [8]. We identify and account for three main factors driving asynchrony in gene expression profiles during a microarray experiment, namely: (1) diversity of infection time; (2) diversity of growth rate; and (3) the emergence of early stage parasites intermixed with later stage parasites. By including these considerations into our computational framework, we are able to more accurately reconstruct single cell, global gene expression profiles, from which the temporal expression profile for individual protein kinases are determined.

Methods

A publicly available microarray data set on three P. falciparum strains (HB3, 3D7 and Dd2) is used in this work [4, 13]. In the HB3 data set, the expression profiles of 4345 oligonucleotide sequences have been measured at 48 time points with 1 hour interval. The 23^rd and 29^th data points are missing for all oligonucleotide sequences. The oligonucleotide sequences associated with protein kinases involved in the P. falciparum life cycle are retrieved from the PlasmoDB database. Some protein kinases have several unique oligonucleotide sequences. For those protein kinases, an average trace is calculated from the curves associated with each oligonucleotide sequence. In this fashion, the gene expression profiles of 65 protein kinases are collected from the data set of HB3. Along the same lines, 52 protein kinases and 51 protein kinases are collected from the data set 3D7 and Dd2 respectively. In the rest of this section, we present a computational method to extract the intrinsic gene expression pattern from the microarray data.

Gene expression levels obtained in microarray experiments are aggregates across many individual iRBCs. First, we will derive a set of linear equations (15) that relates the microarray data to the intrinsic expression pattern of individual iRBC. Next, by solving the corresponding linear inverse problem (16), we reconstruct the expression pattern. All symbols used in this paper are explained in the Table 1.

Table 1 Explanations of symbols

Analysis of decaying synchrony

We assume that three main factors drive the iRBCs out of synchronization in the microarray experiment, namely: (A) diversity of infection time, (B) diversity of growth rate, and (C) the emergence of early stage parasites intermixed with later stage parasites. These are discussed below.

Diversity of infection time

As illustrated in Figure 1, the invasion of RBCs does not occur simultaneously. In the microarray experiment [4], late-stage schizonts are synchronized by six sorbitol treatments on three generations. Prior to the first microarray time point, fresh RBCs are infected by late-stage schizonts within two hours, raising the parasitemia from 5% to 16%. After the invasion period, 80% of parasites are in the ring stage. Let M stand for the total number of cells in the media. In other words, 5% × M of schizonts infect 80% × 16% × M of ring during two-hour invasion, remaining 20% × 16% × M of schizonts are still alive after the invasion period. Although the concentration of RBCs is reduced from 14% to 3.3% immediately after the invasion period, RBCs can still be infected as long as schizonts remain. Therefore, a large amount of RBCs will be infected after the two-hour invasion.

Figure 1

Diversity of infection time among iRBCs. The invasion of RBCs does not occur simultaneously. RBCs will be continuously infected as long as schizonts remain.

Let R(t) denote the number of fresh RBCs infected by schizont at time t (hours). In the perfectly synchronized case, R(t) should be a Dirac delta function, which means all iRBCs are simultaneously infected at the same time. In microarray experiment, however, R(t) has a high value during the invasion period, and it maintains positive as long as schizonts remain.

Let S(t) stand for the number of schizonts which infect RBCs at time t. The parameters a _in and a _af denote the average number of RBCs infected by one schizont during and after the infection period respectively. Therefore, the expression of R(t) can be written as:

$$R\left(t\right)=\left\{\begin{array}{c}{a}_{\mathsf{\text{in}}}S\left(t\right),\mathsf{\text{if}}\mathsf{\text{t}}\in \mathsf{\text{[infection period],}}\\ a_{\mathsf{\text{af}}}S\left(t\right),\mathsf{\text{if}}\mathsf{\text{t}}\in \mathsf{\text{[after infection period]}}\cdot \end{array}\right.$$

(1)

At the end of this section, we explain how to estimate the parameters a _in, a _af and function S(t).

Diversity of growth rate

The iRBCs grow at different rates [4]. Consequently, as illustrated in Figure 2, synchrony gradually decays, even if all iRBCs are simultaneously infected at the same time. Let the random variable $\tilde{L}$ indicate the normalized life span of iRBC, which is the quotient of individual life span L' and the average life span L:

$$\tilde{L}=\frac{L^{\prime }}{L}.$$

(2)

Figure 2

Diversity of growth rate among iRBCs. The iRBCs grow at different rates.

$\tilde{L}$ is assumed to follow a normal distribution: $\tilde{L}~N\left(1,{\sigma }^2\right).$ Hence the probability density function $p_{\tilde{L}}\left(l\right)$of normalized life span $\tilde{L}$ can be written as:

$$p_{\tilde{L}}\left(l\right)=\frac{1}{\sqrt{2\pi {\sigma }^2}}{e}^{-\frac{{\left(l-1\right)}^2}{2{\sigma }^2}},$$

(3)

where the value of σ is adjusted to fit the experimental observations. The details will be discussed later in this paper.

Let {f _i (ℓ), ℓ ∈ [0, L]} be the intrinsic gene expression pattern of protein i on one complete life span, where ℓ denotes the cell age of iRBC (hours). Since f _i (ℓ) represents the common pattern shared by individual RBC, the expression profile of individual iRBC is assumed to be $\left\{f_i\left({\ell }_{\mathsf{\text{re}}}\right),{\ell }_{\mathsf{\text{re}}}=\ell /\tilde{L},\ell \in \left[0,L^{\prime }\right]\right\}$, where ℓ_re denotes the rescaled cell age according to its normalized life span $\tilde{L}$ For instance, one iRBC has been infected for ℓ hours, its current expression level on the protein i is written as $f_i\left(\ell /\tilde{L}\right),$ where $\tilde{L}$ is the normalized life span of the corresponding RBC.

The emergence of early stage parasites intermixed with later stage parasites

Due to the diversity of growth rate, a few iRBCs can reach the late stage of schizont early. As a result, those fast-growing iRBCs can infect additional fresh RBCs, as illustrated in Figure 3. This phenomenon has been observed in experiments [4]. Let S _f (t) denote the number of fast-growing iRBCs which reach end of their life span at time t. Let R _f (t) be the number of RBCs which are infected by these iRBCs at time t. We assume that R _f (t) is proportional to S _f (t):

$$R_f\left(t\right)=a_{\mathsf{\text{af}}}{S}_f\left(t\right).$$

(4)

Figure 3

The emergence of early stage parasites intermixed with later stage parasites. Additional fresh RBCs will be infected by fast-growing iRBCs which reach the late stage of schizont early.

The invasion factor a _af stands for the average number of fresh RBCs that will be infected by one schizont after the invasion period.

As shown in (2), the normalized life span $\tilde{L}$ stands for the quotient of individual life span L' and the average life span L. The number of iRBCs that start and end their life span at time t are denoted as R(t) and S _f (t) respectively. Given the probability density function $p_{\tilde{L}}\left(l\right)$ of normalized life span $\tilde{L}$, the number of iRBCs that have reached the end of their life span at time t can be written as follow:

$$\begin{array}{ll}\hfill {\int }_{-\infty }^t{S}_f\left(t^{\prime }\right)\mathsf{\text{d}}{t}^{\prime }& ={\int }_{-\infty }^{+\infty }R\left(t^{\prime }\right)P_{\tilde{L}}\left(\tilde{L}<\frac{t-t^{\prime }}{L}\right)\mathsf{\text{d}}{t}^{\prime }\\ ={\int }_{-\infty }^{+\infty }R\left(t^{\prime }\right){\int }_{-\infty }^{\frac{t-t^{\prime }}{L}}{p}_{\tilde{L}}\left(l\right)\mathsf{\text{d}}l\mathsf{\text{d}}{t}^{\prime }.\end{array}$$

(5)

Therefore, the expression of S _f (t) can be derived from (5) as:

$$\begin{array}{ll}\hfill S_f\left(t\right)& =\frac{\mathsf{\text{d}}}{\mathsf{\text{d}}t}{\int }_{-\infty }^t{S}_f\left(t^{\prime }\right)\mathsf{\text{d}}{t}^{\prime }\\ =\frac{1}{L}{\int }_{-\infty }^{+\infty }R\left(t^{\prime }\right)p_{\tilde{L}}\left(\frac{t-t^{\prime }}{L}\right)\mathsf{\text{d}}{t}^{\prime }.\end{array}$$

(6)

Therefore, we have the expression of R _f (t) by substituting (6) into (4):

$$R_f\left(t\right)=\frac{a_{\mathsf{\text{af}}}}{L}{\int }_{+\infty }^{+\infty }R\left(t\prime \right)p_{\tilde{L}}\left(\frac{t-t\prime }{L}\right)\mathsf{\text{d}}t\prime .$$

(7)

Simulation of iRBCs population distribution

Let N(t) denote the total number of iRBCs at time t. N(t) consists of 3 parts: the late-schizonts that infect fresh RBCs during the infection period and have not yet burst at time t, the first generation of iRBCs (infected by late-schizonts around infection period) that have not yet reached the end of their life span at time t, and second generation of iRBCs (infected by fast-growing iRBCs) that have not yet reached the end of their life span at time t. Therefore, as illustrated in Figure 4, we can decompose N(t) as:

$$N\left(t\right)={\int }_t^{+\infty }S\left(t^{\prime }\right)\mathsf{\text{d}}{t}^{\prime }+{\int }_{+\infty }^{t}R\left(t^{\prime }\right)p_{\tilde{L}}\left(\tilde{L}>\frac{t-t^{\prime }}{L}\right)\mathsf{\text{d}}{t}^{\prime }+{\int }_{+\infty }^t{R}_f\left(t^{\prime }\right)p_{\tilde{L}}\left(\tilde{L}>\frac{t-t^{\prime }}{L}\right)\mathsf{\text{d}}{t}^{\prime },$$

(8)

Figure 4

The total number of iRBCs N ( t ) at time t. (A) The late-schizonts that infect fresh RBCs during infection period and have not yet burst at time t. (B) The first generation of iRBCs (infected by late-schizonts around infection period) that have not yet reached the end of their life span at time t. (C) Second generation of iRBCs (infected by fast-growing iRBCs) that have not yet reached the end of their life span at time t.

where S(t) stands for the number of late-schizonts bursts at time t, R(t) denotes the number of RBCs infected by late-schizonts at time t, and R _f (t) is the number of RBCs which are infected by fast-growing iRBCs at time t. The expressions of R(t) and R _f (t) are given by (1) and (7) respectively. The expression of S(t) will be derived in (25).

Let N(t, ℓ_re) stand for the number of iRBCs that reach the rescaled cell age ℓ_re at time t. In other words, given the time t, N(t, ℓ_re) indicates the distribution of iRBCs over a complete life span of iRBC (see Figure 5). Therefore, N(t, ℓ_re) and N(t) satisfy following relation:

$$N\left(t\right)={\int }_{-\infty }^{L}N\left(t,{{\ell }^{\prime }}_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{{\ell }^{\prime }}_{\mathsf{\text{re}}}.$$

(9)

Figure 5

iRBCs population distributions change over time. N(t, ℓ_re) with σ = 0.1. It indicates the number of iRBCs that reach the rescaled cell age ℓ_re at time t.

Hence, as illustrated in Figure 6, the number of iRBCs that have not yet reached the rescaled cell age ℓ_re time t can be expanded from (8) as follows:

$$\begin{array}{ccc}\hfill {\int }_{-\infty }^{{\ell }_{\mathsf{\text{re}}}}N\left(t,{{\ell }^{\prime }}_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{{\ell }^{\prime }}_{\mathsf{\text{re}}}\hfill & \hfill =\hfill & \hfill {\int }_t^{+\infty }S\left(t^{\prime }\right)P_{\tilde{L}}\left(\tilde{L}<\frac{t^{\prime }-t}{L-{\ell }_{\mathsf{\text{re}}}}\right)\mathsf{\text{d}}{t}^{\prime }+{\int }_{-\infty }^{t}R\left(t^{\prime }\right)P_{\tilde{L}}\left(\tilde{L}>\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}\right)\mathsf{\text{d}}{t}^{\prime }\hfill \\ \hfill +\hfill & {\int }_{-\infty }^t{R}_f\left(t^{\prime }\right)P_{\tilde{L}}\left(\tilde{L}>\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}\right)\mathsf{\text{d}}{t}^{\prime }\hfill \\ \hfill =\hfill & {\int }_t^{+\infty }S\left(t^{\prime }\right){\int }_{-\infty }^{\frac{t^{\prime }-t}{L-{\ell }_{\mathsf{\text{re}}}}}{p}_{\tilde{L}}\left(\ell \right)\mathsf{\text{d}}l\mathsf{\text{d}}{t}^{\prime }+{\int }_{-\infty }^{t}R\left(t^{\prime }\right){\int }_{\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}}^{+\infty }{p}_{\tilde{L}}\left(l\right)\mathsf{\text{d}}l\mathsf{\text{d}}{t}^{\prime }\hfill \\ \hfill +\hfill & {\int }_{-\infty }^t{R}_f\left(t^{\prime }\right){\int }_{\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}}^{+\infty }{p}_{\tilde{L}}\left(l\right)\mathsf{\text{d}}l\mathsf{\text{d}}{t}^{\prime }.\hfill \end{array}$$

(10)

Figure 6

The number of iRBCs that have not reached the rescaled cell age ℓ _re at time t. (A) The late-schizonts that infect fresh RBCs around infection period and have not reached the rescaled cell age ℓ_re at time t. (B) The first generation of iRBCs that have not reached the rescaled cell age ℓ_re at time t. (C) Second generation of iRBCs that have not reached the rescaled cell age ℓ_re at time t.

Therefore, N(t, ℓ_re) can be derived from (10) as follows:

$$\begin{array}{ll}\hfill N\left(t,{\ell }_{\mathsf{\text{re}}}\right)& =\frac{\mathsf{\text{d}}}{\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}}{\int }_{-\infty }^{{\ell }_{\mathsf{\text{re}}}}N\left(t,{{\ell }^{\prime }}_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}^{\prime }\\ ={\int }_t^{+\infty }S\left(t^{\prime }\right)p_{\tilde{L}}\left(\frac{t^{\prime }-t}{L-{\ell }_{\mathsf{\text{re}}}}\right)\frac{t^{\prime }-t}{{\left(L-{\ell }_{\mathsf{\text{re}}}\right)}^2}\mathsf{\text{d}}{t}^{\prime }+{\int }_{-\infty }^{t}R\left(t^{\prime }\right)p_{\tilde{L}}\left(\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}\right)\frac{t-t^{\prime }}{{{\ell }_{\mathsf{\text{re}}}}^2}\mathsf{\text{d}}{t}^{\prime }\\ +{\int }_{-\infty }^t{R}_f\left(t^{\prime }\right)p_{\tilde{L}}\left(\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}\right)\frac{t-t^{\prime }}{{{\ell }_{\mathsf{\text{re}}}}^2}\mathsf{\text{d}}{t}^{\prime }.\end{array}$$

(11)

By substituting the expression of R _f (t) (7) into (11), N(t, ℓ_re) can be written as:

$$\begin{array}{ll}\hfill N\left(t,{\ell }_{\mathsf{\text{re}}}\right)& ={\int }_t^{+\infty }S\left(t^{\prime }\right)p_{\tilde{L}}\left(\frac{t^{\prime }-t}{L-{\ell }_{\mathsf{\text{re}}}}\right)\frac{t^{\prime }-t}{{\left(L-{\ell }_{\mathsf{\text{re}}}\right)}^2}\mathsf{\text{d}}{t}^{\prime }+{\int }_{-\infty }^{t}R\left(t^{\prime }\right)p_{\tilde{L}}\left(\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}\right)\frac{t-t^{\prime }}{{{\ell }_{\mathsf{\text{re}}}}^2}\mathsf{\text{d}}{t}^{\prime }\\ +\frac{a_{\mathsf{\text{af}}}}{L}{\int }_{-\infty }^t{p}_{\tilde{L}}\left(\frac{t-t^{\prime }}{{\ell }_{\mathsf{\text{re}}}}\right)\frac{t-t^{\prime }}{{{\ell }_{\mathsf{\text{re}}}}^2}\mathsf{\text{d}}{t}^{\prime }{\int }_{-\infty }^{+\infty }R\left(t^{″}\right)p_{\tilde{L}}\left(\frac{t^{\prime }-t^{″}}{L}\right)\mathsf{\text{d}}{t}^{″}.\end{array}$$

(12)

Modeling of gene expression level

The gene expression levels obtained in microarray experiments are aggregates across many individual iRBCs. This superposition across iRBCs is modeled by means of a linear system. Let e _i (t) denote the observed expression level of protein i at time t. As discussed earlier, N(t, ℓ_re) denotes the distribution of iRBCs on rescaled cell age ℓ_re and $\left\{f_i\left({\ell }_{\mathsf{\text{re}}}\right),{\ell }_{\mathsf{\text{re}}}=\ell /\tilde{L},\ell \in \left[0.L^{\prime }\right]\right\}$ stands for the gene expression level of individual iRBC according to its rescaled cell age ℓ_re. Therefore, the observed expression level e _i (t) can be written as a integral over one complete life span of iRBCs as follows:

$$e_i\left(t\right)={\int }_0^{L}N\left(t,{\ell }_{\mathsf{\text{re}}}\right)f_i\left({\ell }_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}.$$

(13)

We represent the continuous function f _i (s) and N(t, ℓ_re) as a series of discrete points {f _i (1), f _i (2),..., f _i (T)}, and {N(t,1), N(t,2),...,N(t,T)} respectively. Therefore, e _i (t) can be approximated as:

$$e_i\left(t\right)\approx \sum _{{\ell }_{\mathsf{\text{re}}}=1}^{L}N\left(t,{\ell }_{\mathsf{\text{re}}}\right)f_i\left({\ell }_{\mathsf{\text{re}}}\right)\Delta {\ell }_{\mathsf{\text{re}}}.$$

(14)

In microarray experiments, gene expression levels are measured at a series of discrete time points. Hence, the resulting observed expression level e _i (t) is a series of discrete value. In the dataset HB3, for example, the expression levels are measured every hour over a period of 48 hours. Therefore, a linear system can be derived based on (14):

$$\underset{A}{\underbrace{\left(\begin{array}{c}N\left(1,1\right)N\left(1,2\right)\dots N\left(1,L\right)\\ N\left(2,1\right)N\left(2,2\right)\dots N\left(2,L\right)\\ ⋮⋮\ddots ⋮\end{array}\right)}}\underset{x}{\underbrace{\left(\begin{array}{c}{f}_i\left(1\right)\\ f_i\left(2\right)\\ ⋮\\ f_i\left(L\right)\end{array}\right)}}=\underset{b}{\underbrace{\left(\begin{array}{c}{e}_i\left(1\right)\\ e_i\left(2\right)\\ ⋮\end{array}\right)}}.$$

(15)

The observation matrix A can be calculated by means of the equation (12). The element of matrix A at row t and column ℓ_re denotes the number of iRBCs that reach the rescaled cell age l _re at time point t. The constant vector b stands for the gene expression levels observed in microarray experiment. The unknown variable vector x is the intrinsic gene expression pattern of individual cell. We can find x by solving the discrete linear inverse problem (15).

Reconstruction of synchronized expression profile

For each protein i, the observed expression level e _i is modeled as the superposition of the expression level f _i of individual iRBCs, as described in (15). To solve the described linear inverse problem (15), we minimize an objective function. The objective function contains the squared error (Ax - b)(Ax - b) ^⊤. The intrinsic gene expression pattern of iRBCs is assumed to be a smooth curve. Therefore, we also include square difference as ${\sum }_{k=1}^{L-1}{\left(x_k-x_{k+1}\right)}^2+{\left(x_L-x_1\right)}^2$ in the objection function. We also need to impose the constraint x ≥ 0 because expression levels are positive. In summary, we compute the intrinsic gene expression pattern x as follows:

$$x=\underset{x\ge 0}{\text{arg}\text{min}}\left\{\underset{\mathsf{\text{square error}}}{\underbrace{\left(Ax-b\right){\left(Ax-b\right)}^{\top }}}+c\underset{\mathsf{\text{gradients}}}{\underbrace{\left[\sum _{k=1}^{L-1}{\left(x_k-x_{k+1}\right)}^2+{\left(x_L-x_1\right)}^2\right]}}\right\}.$$

(16)

We solve (16) numerically by means of the function f mincon in the Optimization Toolbox of Matlab (MATLAB 7.9, The MathWorks Inc.)

Prediction of protein kinases

As we discussed earlier, an increased gene expression level before a cell stage transition is regarded as a sign that the corresponding protein kinase is involved in that stage transition. We denote the likelihood that protein kinases i is involved in regulating the stage transition s by L _i (s). Let T _s stand for the time point when stage transition s occur. The expression of likelihood L _i (s) is derived in following.

Let ${\tilde{f}}_i\left(\ell \right)$ be the normalized expression pattern such that its integral over one iRBC life span is equal to 1:

$${\tilde{f}}_i(\ell )=\frac{f_i(\ell )}{{\displaystyle {\sum }_{\ell =1}^L{f}_i(\ell )\Delta l}},$$

(17)

where ℓ = 1,..., L. Since protein kinases are often regulated at translational and post-translational levels [5], L _i (s) is estimated as the maximum sum of W consecutive points of ${\tilde{f}}_i\left(\ell \right)$ that appear H hours prior to T _s .

$$L_i\left(s\right)=\underset{T_s-H\le {\ell }_1\le T_s-W+1}{\text{max}}\left\{\sum _{\ell ={\ell }_1}^{{\ell }_1+W-1}{\tilde{f}}_i\left(\ell \right)\right\}.$$

(18)

Consequently, protein kinase i can be prioritized by its likelihood L _i (s) of being involved in a stage transition.

Parameter estimation

The proposed method for reconstructing expression patterns consists of two main steps. First, the linear system described in (15) is built by calculating the cell age distribution N(t, ℓ_re) (12) for each element of the observation matrix A. Second, the expression pattern is reconstructed by solving the discrete linear inverse problem (16). Therefore, the performance of our method mainly depends on the calculation of N(t, ℓ_re). As shown in equation (12), N(t, ℓ_re) is dominated by three sets of parameters: the infection factors a _in and a _af, the burst rate of schizonts S(t), and the standard deviation σ of normalized life span $\tilde{L}$.

The infection factors a _in, a _af, and burst rate of schizonts S(t) can be accurately estimated from parasitemia and percent representation of iRBC observed at each time point. As shown in Figure 7, the percent representation of iRBC is available at each time point [4]. However, the parasitemia is given only at two time points, one during and one after the infection period. In this section, we explain how we estimate all parameters from the specification of the microarray experiments [4].

Figure 7

Percent representation of iRBCs observed in experiment. The boundaries (T ₁, T ₂ and L) between the three stages (ring, trophozoite, and schizont) are approximated from the percent representation of iRBCs observed at every time point in experiment [4].

Infection factors

In the proposed model, the infection factors, a _in and a _af denote the number of RBC that can be infected by one bursted schizont during the infection period and after the period respectively. As discussed earlier, 5% × M schizonts infect 80% ×16% × M of ring cell during invasion period, and the remaining 20% × 16% × M schizonts are still alive after the invasion period. Therefore, the value of a _in can be deduced as follows:

$$a_{\mathsf{\text{in}}}=\frac{16\%\times 80\%\times M}{5\%\times M-16\%\times 20\%\times M}\approx 7.11.$$

(19)

After the invasion period, the cell concentration (both fresh and infected) is reduced from 14% to 3.3%. To estimate the value of a _af, we propose a simple model to describe how infection factor is influenced by the cell concentration.

At the end of the life span, the cell membrane of schizont bursts, and merozoites are released to infect other RBCs [14]. We denote by p the probability that a merozoite does not infect any RBC, and let n represent the average number of merozoites released by one schizont. The value of n is estimated to be 14, because one schizont usually contains 12 to 16 merozoites [14]. Therefore, the average number of RBCs infected by one bursted schizont is given by n(1 - p).

In the following, we derive an expression for the probability p. First, we assume that each merozoite can travel a certain space after it is released from schizont, and let V be the volume of this space. Second, we also assume that if a RBC appears in the space indicated by V, it will be immediately infected by the corresponding merozoite. Let V _total denote the volume of whole media. Let m stand for the total number of RBCs remaining in the media. Hence, the value of p can be estimated as the probability that none of the m RBCs appear in the space V. If the RBCs are uniformly distributed in V _total, it follows:

$$p={\left(\frac{V_{\mathsf{\text{total}}}-V}{V_{\mathsf{\text{total}}}}\right)}^m.$$

(20)

The concentration of RBCs is reduced by adjusting the volume of culture V _total from 1000 milliliter to 4500 milliliter after the infection period [4]. Therefore, the expression of a _in and a _af can be obtained by substituting (20) into n(1 - p) as follows:

$$\left\{\begin{array}{c}{a}_{\mathsf{\text{in}}}=14\left(1-{\left(\frac{1000-V}{1000}\right)}^m\right)\\ a_{\mathsf{\text{af}}}=14\left(1-{\left(\frac{4500-V}{4500}\right)}^{m^{\prime }}\right).\end{array}\right.$$

(21)

At the beginning of the experiment, the culture is initialized by 115.0 milliliter of purified RBC [4]. Human blood has 4 to 6 million RBC per microliter (cubic millimeter), and the corresponding hematocrit is about 45% [15]. Hence, the number of RBC m before the infection period can be estimated as:

$$m=\frac{5,000,000\times 115,000}{45\%}\approx 1.28\times 10^{12}.$$

(22)

After the infection period, the parasitemia is raised from 5% to 16%. Therefore, the number of RBC before the infection period m and after the infection period m' satisfy following relation:

$$\frac{m}{100\%-5\%}=\frac{m^{\prime }}{100\%-16\%}.$$

(23)

Consequently, the number of RBC after the infection period m' can be estimated. By substituting the value of a _in, m and m' into (21), the value of V and a _af can be deduced as follows:

$$\left\{\begin{array}{c}V\approx 5.55\times 10^{-10}\mathsf{\text{millilitre,}}\\ a_{\mathsf{\text{af}}}\approx 1.82.\end{array}\right.$$

(24)

The derivation of V is related to the concept of mean free paths in physics, roughly the diffusion length for the first binding event or the lifetime [16].

Burst rate of schizonts

The number of schizonts infecting RBCs at time t is denoted by S(t). As discussed earlier, 5% × M - 16% × 20% × M schizonts burst in the two-hour infection period, leaving 16% × 20% × M schizonts alive till around 7 hours after the invasion period (8th microarray data point). In other words, 36% of schizonts burst in the first two hours followed by the remaining 64% of schizonts which burst within the next 7 hours. Therefore, S(t) is approximated as a piecewise linear function whose integral on these two periods is 36 and 64 respectively. The value of S(t) in the first two hours is also assumed to be a constant. As shown in Figure 8, consequently, the expression of S(t) can be written as follows:

$$s\left(t\right)=\left\{\begin{array}{c}\hfill 18,\hfill \\ \hfill -2.54t+20.54,\hfill \\ \hfill 0,\hfill \end{array}\begin{array}{c}\hfill \mathsf{\text{if}}-1\le t<1,\hfill \\ \hfill \mathsf{\text{if}}1\le t\le 8.1,\hfill \\ \hfill \mathsf{\text{otherwise}}\mathsf{\text{.}}\hfill \end{array}\right.$$

(25)

Figure 8

Burst rate of schizonts. The number of schizonts S(t) infecting RBCs at time t is approximated as a piecewise linear function according to experimental observations.

Standard deviation of normalized life span

Given a _in, a _af, and S(t), the iRBCs population distributions N(t, ℓ_re) (12) depend on the probability density function $p_{\tilde{L}}$(l) of normalized life span $\tilde{L}$ (3). The normalized life span $\tilde{L}$ has been assumed to follow a normal distribution. By definition, the mean of $\tilde{L}$ is equal to 1. The standard deviation σ is estimated such that the resulting iRBCs population distributions N(t, ℓ_re) are in agreement with the experimental data.

Over one complete life cycle, iRBCs go through three life stages, i.e., ring, trophozoite and schizont, as follows:

$$\mathsf{\text{iRBC}}\mathsf{\text{stage}}=\left\{\begin{array}{c}\hfill \mathsf{\text{ring,}}\hfill \\ \hfill \mathsf{\text{trophozoite,}}\hfill \\ \hfill \mathsf{\text{schizont,}}\hfill \end{array}\begin{array}{c}\hfill \mathsf{\text{if}}\mathsf{\text{0}}\le {\ell }_{\mathsf{\text{re}}}<T_1,\hfill \\ \hfill \mathsf{\text{if}}{T}_1\le {\ell }_{\mathsf{\text{re}}}<T_2,\hfill \\ \hfill \mathsf{\text{if}}{T}_2\le {\ell }_{\mathsf{\text{re}}}\le L.\hfill \end{array}\right.$$

(26)

As illustrated in Figure 7, the value of T ₁, T ₂ and L are approximated from the percent representation of iRBCs. Given the boundaries between the three stages, iRBCs population distributions N(t, ℓ_re) can represent the percent representation of iRBCs at each time point as follows. By definition, N(t, ℓ_re) stands for the number of iRBCs that reach the rescaled cell age l _re at time t. Hence, the number of iRBC at a specific stage (ring, trophozoite, and schizont) can be estimated as a integral of N(t, ℓ_re) on the rescaled cell age l _re . Therefore, the percent of ring cell at time t can be calculated as follows:

$$\mathsf{\text{percent of ring at time }}t\mathsf{\text{ = }}\frac{\stackrel{\mathsf{\text{ring}}}{\stackrel{⏟}{{\int }_0^{T_1}N\left(t,{\ell }_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}}}}{\underset{\mathsf{\text{ring}}}{\underbrace{{\int }_0^{T_1}N\left(t,{\ell }_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}+}}\underset{\mathsf{\text{trophozoite}}}{\underbrace{{\int }_{T_1}^{T_2}N\left(t,{\ell }_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}+}}\underset{\mathsf{\text{schizont}}}{\underbrace{{\int }_{T_2}^{L}N\left(t,{\ell }_{\mathsf{\text{re}}}\right)\mathsf{\text{d}}{\ell }_{\mathsf{\text{re}}}.}}}$$

(27)

The percent of trophozoite and schizont can be calculated similarly.

The distribution N(t, ℓ_re) for σ = 0.1 is shown in Figure 5. The percent of ring, trophozoite, and schizont calculated from N(t, ℓ_re) is illustrated in Figure 9. The corresponding calculated percent representation of iRBCs depends on the value of σ. Consequently, we tune the value of σ such that percent representation observed in experiment, as illustrated in Figure 7, coincides with the simulated percent representation, as illustrated in Figure 9.

Figure 9

Percent representation of iRBCs calculated from iRBCs population distributions. The percent of iRBCs (ring trophozoite, and schizont) calculated from iRBCs population N(t, ℓ_re) with σ = 0.1.

Results on synthetic data

In this section, the proposed method is evaluated on synthetic microarray data.

Generate synthetic microarray data

The data obtained in microarray experiments are aggregates across many cells. This superposition across cells is modeled by means of a linear system (15). The matrix A in that linear system depends on the parameters (a _in, a _af, σ, and S(t)), estimated from microarray experiments [4]. Synthetic microarray data is generated by substituting the expression pattern of individual iRBC x into the linear system. As shown in Figure 10, we generate synthetic microarray data for four expression patterns (A, B, C and D). Each of them simulates a protein kinase that serves a specific function in cell stage transition, and hence has high expression level in a short time period before a cell stage transition.

Figure 10

Synthetic microarray data generated from known expression patterns. Synthetic microarray data are generated for four expression patterns (A, B, C and D).

The microarray experiments measure the superposition across cells. Due to the decay of synchrony among cells, the intrinsic expression pattern is blurred in microarray data. By comparing the intrinsic expression patterns and synthetic microarray data shown in Figure 10, we demonstrate how decaying synchrony can blur the expression pattern.

Tolerance to signal noise and missing data points

Due to the assay complexity in microarray experiments, signal noise [17] and missing time points [18] are common phenomena observed in real microarray data. To evaluate the performance of our method under similar conditions, expression patterns are reconstructed from synthetic microarray data contaminated by Gaussian white noise and/or missing data points. The resulting expression patterns are compared to the original expression patterns (A, B, C and D).

As shown in Figure 11, the synthetic microarray data of expression pattern (A) are contaminated by Gaussian white noise. The signal to noise ratios (SNR) of the resulting microarray data are 10, 15 and 20 dB respectively. To calculate matrix A, we choose the same parameter values (a _in, a _af, σ, and S(t)) as used to generate synthetic microarray data. In the next section, we will investigate how the results vary with other choice of those parameter values. The resulting matrix A is substituted into the linear system described in equation (15), and the expression pattern is reconstructed by solving the corresponding linear inverse problem (16). The resulting expression patterns shown in Figure 11 are quite similar to the original expression pattern (A) given in Figure 10.

Figure 11

Microarray data contaminated by Gaussian white noise. Gaussian white noise is used to contaminate the synthetic microarray data of expression pattern (A). The reconstructed patterns are quite similar to the original expression pattern (A).

Along the same lines, the tolerance to missing data points is evaluated using synthetic microarray data with missing data points. We randomly remove data points (10, 20 and 30) from the microarray data of expression pattern respectively. Then expression pattern are reconstructed from contaminated microarray data. As shown in Figure 12, the proposed method generates a reliable estimate for the original expression pattern (A), even when 30 of 48 time points are missing in the microarray data.

Figure 12

Microarray data contaminated by random missing data points. Data points (10, 20 and 30) are randomly remove from the microarray data of expression pattern. The expression pattern can be reliably reconstructed, even when 30 of 48 time points are missing in the microarray data.

Real microarray data are usually contaminated by both signal noise and missing data points. Therefore, as shown in Figure 13, we conduct additional experiments on synthetic microarray data contaminated by both Gaussian white noise (SNR = 10, 20, and 30 dB) and 20 missing data points. The results of Figure 13 suggest that the proposed method can reliably reconstruct the expression pattern when the microarray data is simultaneously contaminated by both signal noise and missing data points.

Figure 13

Microarray data contaminated by both Gaussian white noise and random missing data points. The proposed method can reliably reconstruct the expression pattern when the microarray data is simultaneously contaminated by both signal noise (SNR = 10, 20, and 30 dB) and missing data points (20).

Sensitivity to model parameters

The performance of our method mainly depends on the calculation of cell age distribution N(t, ℓ_re). As shown in equation (12), N(t, ℓ_re) depends on three sets of parameters: the infection factors a _in and a _af, the bufirst rate of schizonts S(t), and the standard deviation σ of normalized life span $\tilde{L}$. In the experiments of the previous section, the same parameters are used to generate synthetic data and to reconstruct expression patterns. However, in real microarray data, these parameters are unknown, and need to be estimated from experimental observations. The infection factors a _in, a _af, and burst rate S(t) of schizonts can be accurately calculated from parasitemia and percent representation of iRBC. It is difficult to precisely estimate the standard deviation σ of normalized life span $\tilde{L}$. In this section, we investigate how results vary with regards to the choice of σ.

Synthetic microarray data is generated from know expression pattern with σ = 0.1. The synthetic microarray data is contaminated by Gaussian white noise (SNR = 10 dB) and 20 missing data points. The expression pattern is reconstructed with various values of σ (σ=0.05, 0.1, or 0.15). As shown in Figure 14, the reconstructed expression patterns only slightly change as the value of σ varies. Consequently, the reconstruction is robust to the choice of σ.

Figure 14

Sensitivity to model parameter σ. Synthetic microarray data is generated from know expression pattern with σ = 0.1. The synthetic microarray data is contaminated by Gaussian white noise (SNR = 10 dB) and 20 missing data points. The reconstructed expression patterns only slightly change as the value of σ varies (σ =0.05, 0.1, or 0.15).

Results on real data

Expression patterns are reconstructed for 68 protein kinases collected from the microarray data of P. falciparum (HB3, 3D7 and Dd2) [4, 13]. Reconstructed expression patterns are substituted into the equation (18) to estimate the likelihood of each protein kinase being associated with a specific IDC transition, and hence, contributing to effecting either the stage transition itself or a process(es) needed in the subsequent stage. Data for the three broad transitions analyzed, namely ring-to-trophozoite, trophozoite-to-schizont and schizont-to-ring, are summarized in Additional Files 1, 2 and 3, respectively.

A primary motivation for developing this computational framework is to prioritize gene candidates with potentially important stage-dependent functions for detailed downstream experimental analysis of gene function. In the ring-to-trophozoite analysis, several members of the largely unstudied FIKK protein kinases emerge with relatively high probabilities of mediating important biology during this developmental transition. Several of these protein kinases are targeted to the infected RBC cytosol/membrane [19]. Two [MAL7P1.144, PFL0040c] have been previously studied using gene knockout approaches [20]. While non-essential to blood stage parasite growth and survival, these proteins help mediate the increased rigidity of infected RBCs observed in trophozoite stage parasites. Presumably, this requires modulation of the RBC cytoskeleton through a combination of RBC cytoskeletal and/or exported parasite protein phosphorylation and increased interactions between these [20]. The analysis here suggests that the other highly ranked family members (see Additional File 1) could also be mediating important yet unknown biology at this ring-trophozoite transition.

Interestingly, in the trophozoite-schizont and schizont-ring analyses, a number of protein kinases previously established to be essential and in some cases implicated in directly impacting the transition emerge with the highest probability rankings (Additional Files 2 and 3) [6–8, 12]. In the schizont-ring analysis, for example, PFB0815 has previously been implicated in parasite motility/invasion/egress [6], and PF13_-0211 in parasite egress from the RBC [7]. Furthermore, MAL13P1.278 is an essential Aurora kinase (Pfark1) that associates with spindle pole bodies during parasite schizogeny and is implicated in cell cycle regulation [8]. Overall, the top ranked genes in the trophozoite-schizont and schizont-ring analyses are highly enriched for protein kinases previously established to be essential to parasite survival. Therefore, it will be intriguing to experimentally examine the other highly ranked genes in both transition categories [PFI1415w, PF11_-0079, PF10_-0380, PF11_-0510, PFI0095c, PFE0045c and PFC0945w] as these may also play critical roles in regulating parasite development during these stages and could provide new opportunities for antimalarial drug development.

Conclusions

This study proposes a new methodology to reconstruct intrinsic expression patterns from microarray data. We derive a linear system that relates the microarray data to the expression patterns. By solving the corresponding linear inverse problem, the expression patterns are reconstructed. The experiments conducted on synthetic data suggest that the proposed method can reliably reconstruct the expression pattern, even though both signal noise and missing data points contaminate the microarray data. By applying this method to P. falciparum microarray data, protein kinases are prioritized in terms of their likelihood of being involved in regulating some aspect(s) of the IDC. Indeed, the results of our analyses are supported by previously published experimental data confirming the involvement of several protein kinases in regulating parasite biology at or between developmental stage transitions in the IDC. Our results indicate that experimentally investigating the function of other putative protein kinases not previously studied but ranked highly in our analysis may provide new insights into P. falciparum biology.