Volume 10 Supplement 1

Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2011: Proteome Science

Open Access

Computational synchronization of microarray data with application to Plasmodium falciparum

Proteome Science201210(Suppl 1):S10

https://doi.org/10.1186/1477-5956-10-S1-S10

Published: 21 June 2012

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 [68]. 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 [68]. 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 23rd and 29th 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

Symbols

Explanations

M

the total amount of cells in the media, it consists of RBCs and iRBCs

S(t)

the number of schizonts which infect RBCs at time t

R(t)

the number of fresh RBCs infected by schizont at time t

a in

the average number of RBCs infected by one schizont during the infection period

a af

the average number of RBCs infected by one schizont after the infection period

L ̃

the normalized life span of individual iRBCs

p L ̃ ( l )

the probability density function of normalized life span L ̃

L

the average life span of iRBCs

L'

the individual life span of iRBCs

the cell age of iRBC in hours

{f i (ℓ), ℓ [0, L]}

the intrinsic gene expression pattern of protein i

re

the rescaled cell age according to its normalized life span L ̃

{f i (ℓre), re = L ̃ , ℓ [0, L']}

the gene expression profile of individual iRBC on protein i

S f (t)

the number of fast-growing iRBCs which infect RBCs at time t

R f (t)

the number of RBCs infected by fast-growing iRBCs at time t

N(t)

the total number of iRBCs that have been infected at time t

N(t, ℓre)

the number of iRBCs which reach rescaled cell age ℓre at time t

e i (t)

the observed expression level of protein i at time t

f ̃ i ( )

the normalized gene expression pattern on protein i

s

the transitions between three stages of iRBC: ring, trophozoite, and schizont

L i ( s )

the likelihood that protein kinases i is involved in regulating the stage transition s

T s

the time point when stage transition s occur

n

the average number of merozoites released by one schizont

p

the probability that one merozoites does not infect any RBC

V

the volume of merozoites can travel after it is released from schizont

V total

the volume of whole media

m

the total number of RBCs in the media

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 ( t ) = a in S ( t ) , if t [infection period], a af S ( t ) , if t [after infection period]
(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 L ̃ indicate the normalized life span of iRBC, which is the quotient of individual life span L' and the average life span L:
L ̃ = L L .
(2)
Figure 2

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

L ̃ is assumed to follow a normal distribution: L ̃ ~ N ( 1 , σ 2 ) . Hence the probability density function p L ̃ ( l ) of normalized life span L ̃ can be written as:
p L ̃ ( l ) = 1 2 π σ 2 e - ( l - 1 ) 2 2 σ 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 { f i ( re ) , re = / L ̃ , [ 0 , L ] } , where ℓre denotes the rescaled cell age according to its normalized life span L ̃ For instance, one iRBC has been infected for ℓ hours, its current expression level on the protein i is written as f i ( / L ̃ ) , where 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 ( t ) = a af S f ( t ) .
(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 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 L ̃ ( l ) of normalized life span L ̃ , the number of iRBCs that have reached the end of their life span at time t can be written as follow:
- t S f ( t ) d t = - + R ( t ) P L ̃ L ̃ < t - t L d t = - + R ( t ) - t - t L p L ̃ ( l ) d l d t .
(5)
Therefore, the expression of S f (t) can be derived from (5) as:
S f ( t ) = d d t - t S f ( t ) d t = 1 L - + R ( t ) p L ̃ t - t L d t .
(6)
Therefore, we have the expression of R f (t) by substituting (6) into (4):
R f ( t ) = a af L + + R ( t ) p L ̃ t - t L d t .
(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 ( t ) = t + S ( t ) d t + + t R ( t ) p L ̃ L ̃ > t - t L d t + + t R f ( t ) p L ̃ L ̃ > t - t L d t ,
(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 ( t ) = - L N ( t , re ) d 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:
- re N ( t , re ) d re = t + S ( t ) P L ̃ L ̃ < t - t L - re d t + - t R ( t ) P L ̃ L ̃ > t - t re d t + - t R f ( t ) P L ̃ L ̃ > t - t re d t = t + S ( t ) - t - t L - re p L ̃ ( ) d l d t + - t R ( t ) t - t re + p L ̃ ( l ) d l d t + - t R f ( t ) t - t re + p L ̃ ( l ) d l d t .
(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:
N ( t , re ) = d d re - re N ( t , re ) d re = t + S ( t ) p L ̃ t - t L - re t - t ( L - re ) 2 d t + - t R ( t ) p L ̃ t - t re t - t re 2 d t + - t R f ( t ) p L ̃ t - t re t - t re 2 d t .
(11)
By substituting the expression of R f (t) (7) into (11), N(t, ℓre) can be written as:
N ( t , re ) = t + S ( t ) p L ̃ t - t L - re t - t ( L - re ) 2 d t + - t R ( t ) p L ̃ t - t re t - t re 2 d t + a af L - t p L ̃ t - t re t - t re 2 d t - + R ( t ) p L ̃ t - t L d t .
(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 { f i ( re ) , re = / L ̃ , [ 0 . L ] } 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 ( t ) = 0 L N ( t , re ) f i ( re ) d 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 ( t ) re = 1 L N ( t , re ) f i ( re ) Δ 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):
N ( 1 , 1 ) N ( 1 , 2 ) N ( 1 , L ) N ( 2 , 1 ) N ( 2 , 2 ) N ( 2 , L ) A f i ( 1 ) f i ( 2 ) f i ( L ) x = e i ( 1 ) e i ( 2 ) b .
(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 k = 1 L - 1 ( x k - x k + 1 ) 2 + ( x L - x 1 ) 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 = arg min x 0 ( A x - b ) ( A x - b ) square error + c k = 1 L - 1 ( x k - x k + 1 ) 2 + ( x L - x 1 ) 2 gradients .
(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 f ̃ i ( ) be the normalized expression pattern such that its integral over one iRBC life span is equal to 1:
f ˜ i ( ) = f i ( ) = 1 L f i ( ) Δ 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 f ̃ i ( ) that appear H hours prior to T s .
L i ( s ) = max T s - H 1 T s - W + 1 = 1 1 + W - 1 f ̃ i ( ) .
(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 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 1, T 2 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 in = 16 % × 80 % × M 5 % × M - 16 % × 20 % × M 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 = V total - V V total 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:
a in = 14 1 - 1000 - V 1000 m a af = 14 1 - 4500 - V 4500 m .
(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 = 5 , 000 , 000 × 115 , 000 45 % 1 . 28 × 1 0 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:
m 100 % - 5 % = m 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:
V 5 . 55 × 1 0 - 10 millilitre, a af 1 . 82 .
(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 ( t ) = 18 , - 2 . 54 t + 20 . 54 , 0 , if - 1 t < 1 , if 1 t 8 . 1 , otherwise .
(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 L ̃ (l) of normalized life span L ̃ (3). The normalized life span L ̃ has been assumed to follow a normal distribution. By definition, the mean of 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:
iRBC stage = ring, trophozoite, schizont, if 0 re < T 1 , if T 1 re < T 2 , if T 2 re L .
(26)
As illustrated in Figure 7, the value of T 1, T 2 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:
percent of ring at time  t  =  0 T 1 N ( t , re ) d re ring 0 T 1 N ( t , re ) d re + ring T 1 T 2 N ( t , re ) d re + trophozoite T 2 L N ( t , re ) d re . schizont
(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 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 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) [68, 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.

Declarations

Acknowledgements

This work was supported by the Singapore-MIT Alliance (SMA3) Graduate Fellowship granted to Z.W., and DOD ARO (Grant Number W911NF-09-1-0480).

This article has been published as part of Proteome Science Volume 10 Supplement 1, 2012: Selected articles from the IEEE International Conference on Bioinformatics and Biomedicine 2011: Proteome Science. The full contents of the supplement are available online at http://www.proteomesci.com/supplements/10/S1.

Authors’ Affiliations

(1)
Singapore-MIT Alliance for Research and Technology, Centre for Life Sciences
(2)
School of Electrical and Electronic Engineering, Nanyang Technological University
(3)
Department of Biological Engineering, Massachusetts Institute of Technology
(4)
Department of Chemistry, Massachusetts Institute of Technology

References

  1. World Malaria Report, 2011. Geneva, Switzerland: World Health Organization; 2011. [http://www.who.int/malaria/world_malaria_report_2011/en/]
  2. Trampuz A, Jereb M, Muzlovic I, Prabhu RM: Clinical review: Severe malaria. Crit Care 2003,7(4):315–323. [http://dx.doi.org/10.1186/cc2183] 10.1186/cc2183PubMed CentralPubMedView ArticleGoogle Scholar
  3. Le Roch KG, Zhou Y, Blair PL, Grainger M, Moch JK, Haynes JD, De la Vega P, Holder AA, Batalov S, Carucci DJ, Winzeler EA: Discovery of gene function by expression profiling of the malaria parasite life cycle. Science 2003,301(5639):1503–1508. [http://dx.doi.org/10.1126/science.1087025] 10.1126/science.1087025PubMedView ArticleGoogle Scholar
  4. Bozdech Z, Llinás M, Lee B, Wong ED, Zhu J, DeRisi JL: The Transcriptome of the Intraerythrocytic Developmental Cycle of Plasmodium falciparum. PLoS Biol 2003, 1: e5+. [http://dx.doi.org/10.1371/journal.pbio.0000005]PubMed CentralPubMedView ArticleGoogle Scholar
  5. Le Roch KG, Johnson JR, Florens L, Zhou Y, Santrosyan A, Grainger M, Yan SF, Williamson KC, Holder AA, Carucci DJ, Yates JR, Winzeler EA: Global analysis of transcript and protein levels across the Plasmodium falciparum life cycle. Genome Research 2004,14(11):2308–2318. [http://dx.doi.org/10.1101/gr.2523904] 10.1101/gr.2523904PubMed CentralPubMedView ArticleGoogle Scholar
  6. Kato N, Sakata T, Breton G, Le Roch KG, Nagle A, Andersen C, Bursulaya B, Henson K, Johnson J, Kumar KA, Marr F, Mason D, McNamara C, Plouffe D, Ramachandran V, Spooner M, Tuntland T, Zhou Y, Peters EC, Chatterjee A, Schultz PG, Ward GE, Gray N, Harper J, Winzeler EA: Gene expression signatures and small-molecule compounds link a protein kinase to Plasmodium falciparum motility. Nature Chemical Biology 2008,4(6):347–356. [http://dx.doi.org/10.1038/nchembio.87] 10.1038/nchembio.87PubMedView ArticleGoogle Scholar
  7. Dvorin JD, Martyn DC, Patel SD, Grimley JS, Collins CR, Hopp CS, Bright AT, Westenberger S, Winzeler E, Blackman MJ, Baker DA, Wandless TJ, Duraisingh MT: A Plant-Like Kinase in Plasmodium falciparum Regulates Parasite Egress from Erythrocytes. Science 2010,328(5980):910–912. [http://dx.doi.org/10.1126/science.1188191] 10.1126/science.1188191PubMed CentralPubMedView ArticleGoogle Scholar
  8. Reininger L, Wilkes JM, Bourgade H, Miranda-Saavedra D, Doerig C: An essential Aurora-related kinase transiently associates with spindle pole bodies during Plasmodium falciparum erythrocytic schizogony. Molecular microbiology 2011, 79: 205–221. [http://dx.doi.org/10.1111/j.1365–2958.2010.07442.x] 10.1111/j.1365-2958.2010.07442.xPubMed CentralPubMedView ArticleGoogle Scholar
  9. Ward P, Equinet L, Packer J, Doerig C: Protein kinases of the human malaria parasite Plasmodium falciparum: the kinome of a divergent eukaryote. BMC Genomics 2004, 5: 79. [http://dx.doi.org/10.1186/1471–2164–5-79] 10.1186/1471-2164-5-79PubMed CentralPubMedView ArticleGoogle Scholar
  10. Srinivasan N, Krupa A: A genomic perspective of protein kinases in Plasmodium falciparum. Proteins 2005, 58: 180–189. [http://dx.doi.org/10.1002/prot.20278]PubMedGoogle Scholar
  11. Doerig C, Billker O, Haystead T, Sharma P, Tobin AB, Waters NC: Protein kinases of malaria parasites an update. Trends Parasitol 2008,24(12):570–577. [http://dx.doi.org/10.1016/j.pt.2008.08.007] 10.1016/j.pt.2008.08.007PubMedView ArticleGoogle Scholar
  12. Solyakov L, Halbert J, Alam MM, Semblat JP, Dorin-Semblat D, Reininger L, Bottrill AR, Mistry S, Abdi A, Fennell C, Holland Z, Demarta C, Bouza Y, Sicard A, Nivez MP, Eschenlauer S, Lama T, Thomas DC, Sharma P, Agarwal S, Kern S, Pradel G, Graciotti M, Tobin AB, Doerig C: Global kinomic and phospho-proteomic analyses of the human malaria parasite Plasmodium falciparum. Nature Communications 2011, 2: 565+. [http://dx.doi.org/10.1038/ncomms1558]PubMedView ArticleGoogle Scholar
  13. Llinás M, Bozdech Z, Wong ED, Adai AT, DeRisi JL: Comparative whole genome transcriptome analysis of three Plasmodium falciparum strains. Nucleic acids research 2006,34(4):1166–1173. [http://dx.doi.org/10.1093/nar/gkj517] 10.1093/nar/gkj517PubMed CentralPubMedView ArticleGoogle Scholar
  14. Crewe W, Haddock DRW: Parasites and human disease. Wiley; 1985. [http://www.worldcat.org/isbn/9780471010630]Google Scholar
  15. Dean L: Blood groups and red cell antigens. Bethesda, Md. : National Center for Biotechnology Information; 2005. [http://www.ncbi.nlm.nih.gov/books/NBK2261/]Google Scholar
  16. Silbey RJ, Alberty RA, Bawendi MG: Physical Chemistry. John Wiley & Sons; 2005. [http://www.mathematica-journal.com/issue/v9i3/newpublications/ISBN047121504X.html]Google Scholar
  17. Huber W, von Heydebreck A, Sültmann H, Poustka A, Vingron M: Variance stabilization applied to microarray data calibration and to the quantification of differential expression. Bioinformatics 2002,18(suppl 1):S96-S104. [http://bioinformatics.oxfordjournals.org/content/18/suppl_1/S96.abstract] 10.1093/bioinformatics/18.suppl_1.S96PubMedView ArticleGoogle Scholar
  18. Liew AW, Law NF, Yan H: Missing value imputation for gene expression data: computational techniques to recover missing data from available information. Briefings in bioinformatics 2011,12(5):498–513. [http://dx.doi.org/10.1093/bib/bbq080] 10.1093/bib/bbq080PubMedView ArticleGoogle Scholar
  19. Nunes MC, Goldring DP, Doerig C, Scherf A: A novel protein kinase family in Plasmodium falciparum is differentially transcribed and secreted to various cellular compartments of the host cell. Molecular microbiology 2007,63(2):391–403. [http://dx.doi.org/10.1111/j.1365–2958.2006.05521.x]PubMedView ArticleGoogle Scholar
  20. Nunes MC, Okada M, Scheidig-Benatar C, Cooke BM, Scherf A: Plasmodium falciparum FIKK Kinase Members Target Distinct Components of the Erythrocyte Membrane. PLoS ONE 2010,5(7):e11747+. [http://dx.doi.org/10.1371/journal.pone.0011747]PubMed CentralPubMedView ArticleGoogle Scholar

Copyright

© Zhao et al.; licensee BioMed Central Ltd. 2012

This article is published under license to BioMed Central Ltd. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.