A systems model of vesicle trafficking in Arabidopsis pollen tubes A systems model of vesicle trafficking in Arabidopsis pollen tubes

A systems model that describes vesicle trafﬁcking during pollen tube growth in Arabidopsis ( Arabidopsis thaliana ) was constructed. The model is composed of ordinary differential equations that connect the molecular functions of genes expressed in pollen. The current model requires soluble N -ethylmaleimide-sensitive fusion protein attachment protein receptors (SNAREs) and small GTPases, Arf or Rab, to reasonably predict tube growth as a function of time. Tube growth depends on vesicle trafﬁcking that transports phospholipid and pectin to the tube tip. The vesicle trafﬁcking genes identiﬁed by analyzing publicly available transcriptome data comprised 328 genes. Fourteen of them are up-regulated by the gibberellin signaling pathway during pollen development, which includes the SNARE genes SYP124 and SYP125 and the Rab GTPase gene RABA4D. The model results adequately ﬁt the pollen tube growth of both previously reported wild-type and raba4d knockout lines. Furthermore, the difference of pollen tube growth in syp124 / syp125 single and double mutations was quantitatively predicted based on the model analysis. In general, a systems model approach to vesicle trafﬁcking arguably demonstrated the importance of the functional connections in pollen tube growth and can help guide future research directions.

During pollination, when the pollen-stigma interaction is compatible, pollen cells grow a long tube to deliver sperm to the ovule. This polarized growth largely depends on vesicle trafficking that regulates the flow of macromolecules, such as phospholipids and polysaccharides, to the tip of the pollen tube. Studies investigating vesicle trafficking in pollen tubes of tobacco (Nicotiana tabacum) and lily (Lilium longiflorum) revealed the importance of spatial and temporal changes (Moscatelli et al., 2007;Bove et al., 2008;Zonia and Munnik, 2008) as well as the cytoskeletal and calcium regulations (Cai and Cresti, 2009) of vesicles during the period of polarized growth.
On the other hand, transcriptome analysis revealed genes expressed in Arabidopsis (Arabidopsis thaliana) pollen. Four independent studies reported that pollen of Arabidopsis express approximately 30% of genes encoded in the genome (Honys and Twell, 2004;Pina et al., 2005;Borges et al., 2008;Wang et al., 2008). Gene ontology analysis of these transcriptome data demonstrated that the expressed genes in pollen are enriched with vesicle trafficking machinery (Pina et al., 2005). In addition to these studies, transcriptome studies addressing hormone signaling pathways in pollen de-velopment and pollen tube growth have been conducted in Arabidopsis and rice (Oryza sativa; Achard et al., 2004;Cheng et al., 2004;Kaneko et al., 2004;Tsuji et al., 2006;Allen et al., 2007). GA regulates cell growth and development in plants (Silverstone et al., 1997). Several independent studies showed that loss-of-function mutations in the GA signaling pathway impaired pollen development and tube growth in both Arabidopsis and rice. For instance, ga1-3, a GAdeficient Arabidopsis mutant, exhibits a defect in pollen maturation and also shows retarded vegetative growth (Silverstone et al., 1997). In addition, genetic analyses of the GA signaling pathway, namely DELLA (a negative regulator of GA signaling) and GAMYB (a positive regulator of GA signaling) knockouts (Achard et al., 2004;Cheng et al., 2004), and ectopic expression of microRNA (35S:miR159a suppresses the expression of GAMYB in pollen; Achard et al., 2004;Allen et al., 2007) demonstrated that loss-of-function mutations in the GA signaling pathway impaired pollen development and pollen tube growth in Arabidopsis. Similarly, genetic analyses of the GA signaling pathway in rice, namely the GAMYB loss-of-function mutation (gamyb; Kaneko et al., 2004) and miR159 ectopic expression (POsACT1:miR159; Tsuji et al., 2006), showed that the GA signaling pathway regulates rice pollen development as well. Furthermore, analysis of the KAO missense mutation, one of the GA synthesis proteins, has demonstrated defective pollen tube growth in rice (Chhun et al., 2007). These studies suggest that the GA signaling pathway may regulate pollen tube growth in both Arabidopsis and rice.
Pollen cells are enclosed by a phospholipid bilayer (plasma membrane) and a cell wall (polysaccha-rides). The cells also have the endomembrane system, including organelles such as the endoplasmic reticulum (ER), Golgi apparatus, endosome, and vacuole, which are also enclosed by a phospholipid bilayer. During pollen tube growth, several molecules required for pollen tube growth, such as proteins, phospholipids, and polysaccharides, are produced, stored, or transported in these organelles. At the tube tip, phospholipids are incorporated into the plasma membrane while polysaccharides such as pectin are incorporated into the cell wall (Samaj et al., 2006;Lee and Yang, 2008). In order to transport these molecules (cargos) among the organelles and to the tube tip (compartments hereafter), vesicles bud from a compartment membrane and traffic in the cytosol. When a vesicle arrives at its destination, the vesicle membrane fuses to the destination membrane and releases its cargo into or onto the compartment membrane.
Studies on vesicle trafficking in yeast and animals have demonstrated that the minimum requirement for protein budding and vesicle fusion is small GTPases, namely Arf or Rab, and soluble N-ethylmaleimidesensitive fusion protein attachment protein receptors (SNAREs; Seabra and Wasmeier, 2004;Behnia and Munro, 2005;Hofmann et al., 2006;Cai et al., 2007;Markgraf et al., 2007). During vesicle budding, a small GTPase located on a compartment membrane recruits machinery proteins, such as coat proteins, from the cytosol. These machinery proteins generate vesicles from the compartment membrane. At the same time, the small GTPase initiates the recruitment of vesicle (v)-SNAREs localized in the compartment membrane and cargos into the vesicle through a cascade of protein interactions. While a vesicle travels within the cytosol, a small GTPase on the vesicle membrane may be converted to the other small GTPase that defines the vesicle destination through a series of protein interactions (Grosshans et al., 2006). Meantime, v-SNARE and the cargo remain in the vesicle. When the vesicle reaches its destination, v-SNAREs interact with target (t)-SNAREs in the compartment membrane (Sorensen et al., 2006). This interaction facilitates vesicle fusion using energy generated between v-SNAREs and t-SNAREs (Stein et al., 2009). To date, more than 15 different gene families whose products are part of a complex machinery of vesicle trafficking have been identified (Schmid and McMahon, 2007). Although not included in the vesicle trafficking family in the strict sense, cytoskeletal proteins such as actins, myosins, microtubules and dyneins are also known to play an important role in vesicle transportation within the cytosol in yeast and animals (Ridley, 2006). Because plant genomes encode homologs of all of these genes, it has been suggested that the basic mechanism of vesicle trafficking in plants is similar to that of yeast and animals (Sanderfoot and Raikhel, 2003;Samaj et al., 2005;Lam et al., 2007;Paul and Frigerio, 2007;Otegui and Spitzer, 2008;Robinson et al., 2008;Rojo and Denecke, 2008).
Direct evidence that shows the importance of vesicle trafficking machinery in pollen tube growth has been mainly reported using gene knockout in Arabidopsis, where pollen tube growth is impaired either in vitro or in vivo. These include knockouts of small GTPase (Szumlanski and Nielsen, 2009), phosphatidylinositol (PI) kinase (Ischebeck et al., 2008;Sousa et al., 2008), vesicle tethering factor (Cole et al., 2005), vacuolar membrane protein (Fujiki et al., 2007), and vesicle coat protein (Van Damme et al., 2006). However, the fundamental question of how these gene products lead to pollen tube growth remains unanswered.
Systems biology views a biological process as a complex system rather than trying to reduce or simplify the process into discrete elements (Kitano, 2002). A central goal of systems biology is to construct a mathematical model that predicts the outcome of a biological process of interest. Because genes expressed in Arabidopsis pollen tubes have been identified and the system outcome (tube growth) can be quantitatively measured (both in vitro and in vivo) as a function of time, this process represents an ideal arena in which to take full advantage of a systems biology approach and the rich genetic materials to answer fundamental questions. Here, we present such a mathematical model that connects the molecular functions of genes expressed in Arabidopsis pollen tubes.

Transcriptome Analysis Elucidates Genes Encoding the Machinery of Vesicle Trafficking in Arabidopsis Pollen
To construct a mathematical model that is capable of being evaluated genetically, we sought to define genes of the vesicle trafficking machinery in Arabidopsis pollen tubes. In Arabidopsis, four independent studies on transcriptional profiling of pollen using the same microarray format (Affymetrix Arabidopsis ATH1 Genome Arrays) have been reported (Honys and Twell, 2004;Pina et al., 2005;Borges et al., 2008;Wang et al., 2008). The array contained a 22,591-probe set, which represented approximately 75% of the annotated genes in the Arabidopsis genome (Wang et al., 2008). All studies came to a similar conclusion that approximately 30% of the genes probed on the array (about 6,500 genes) were expressed in mature pollen or pollen tubes.
This categorization revealed that 49% (161 of 328) of the genes are present and 51% (167 of 328) of the genes are absent in pollen (Fig. 1). Of these present genes in pollen, 13% (42 of 328) are enriched or selectively expressed in pollen (Fig. 1). This suggests that a unique set of genes, which are not expressed or are expressed at very low levels in other tissues, is involved in vesicle trafficking of Arabidopsis pollen tube growth.

Gene Expression of Selected Small GTPases and SNAREs Is Controlled by the GA Signaling Pathway in Pollen
Next, we sought to identify mutant Arabidopsis lines in which the vesicle trafficking genes expressed in pollen are down-regulated, so that we might be able to find an upstream gene(s) that regulates the expression of vesicle trafficking genes in pollen. To identify mutant lines, we used the 161 vesicle trafficking genes expressed in pollen ( Fig. 1; Supplemental Table S1) as queries in analysis by Genevestigator (Hruz et al., 2008). This analysis revealed two related mutant lines in which the same cluster of vesicle trafficking genes is clearly down-regulated (less than 0.5-fold) in flower tissues (including pollen) compared with the wild type (Table II). The two related mutant lines identified are ga1-3 and 35S:miR159a. We also found that the genes in the cluster were down-regulated in a penta mutant (Cao et al., 2006), in which gal1-3 and the members of the DELLA gene family are knocked out (Table II). The number of genes in the cluster is 14 ( Fig. 1; Table II). Of these, 10 genes belong to the Rab GTPase gene family (and Rab GDI and Arf GAP) or SNARE gene family (Table II). Three genes belong to the exocyst family, and one gene belongs to the PI Table I. Gene families of vesicle trafficking machinery Gene families that are known to compose the vesicle trafficking machinery in yeast and animals are listed. The family members in Arabidopsis were identified by different research groups (indicated in the "Reference" column). The numbers of each family member are listed in the "Members" column. The numbers of the members whose expression is enriched (expressed at least 1.2-fold higher than other tissues) or selective (Affymetrix detection call is present only in pollen) are listed in the "Pollen Enriched or Selective" column. The analysis was based on the data of Honys and Twell (2004), Pina et al. (2005), and Borges et al. (2008). The numbers of members whose expression is down-regulated (less than 0.5-fold) in flower tissues (including pollen) of both ga1-3 (Cao et al., 2006) and 35S:mir159a (Schwab et al., 2005) mutants are listed in the "GA Down" column. The down-regulated genes were identified by Genevestigator V3 (Hruz et al., 2008)  metabolic enzyme family. This suggests that the function of the selected minimum machinery of vesicle trafficking in pollen tube, SNARE and small GTPase, may be controlled by the GA signaling pathway. Exocyst proteins are a subunit of a tethering complex that functions during vesicle fusion in the plasma membranes of yeast and animals (He and Guo, 2009). In Arabidopsis, a knockout mutant of AtEXP70A1 that is absent in pollen (Supplemental Table S1) shows defects in the polar growth of root hairs and stigmatic papillae (Synek et al., 2006). Knocking out genes of several other subunits in the tethering complex causes defects in pollen germination and pollen tube growth (Hala et al., 2008). Moreover, in an Arabidopsis knockout mutant of a P-type ATPase cation pump gene that is required for normal pollen development, male gametogenesis impaired anthers, expression of the three exocyst genes (AtEXO70H3, AtEXO70C1, and AtEXO70C2) is selectively down-regulated (Jakobsen et al., 2005). These results suggest that exocyst proteins may play an important role in pollen tube growth. Knockout analyses in these genes are required to conclude the finding.
The other gene family identified was the PI metabolic enzyme family, PIP5K10. An Arabidopsis knockout mutant of PIP5K4, which is enriched in pollen (Supplemental Table S1) and mainly localizes in the plasma membrane of a pollen tube tip, impairs pollen tube growth (Sousa et al., 2008). Hence, it is reasonable to assume that expression of an enzyme in the same metabolic pathway controls pollen tube growth. Again, knockout analysis in this gene is required to conclude the finding.
We also analyzed changes of the expression levels in these 14 genes during pollen tube growth based on the data published by Wang et al. (2008). We found that the expression levels of these 14 genes are little changed during pollen tube growth (Wang et al., 2008), suggesting that the vesicle trafficking may be established rapidly before or during pollen germination. Furthermore, we found that 10 of the 14 genes, which are all enriched and selectively expressed in pollen in the data of Pina et al. (2005), are up-regulated Table II. Vesicle trafficking genes that are down-regulated in Arabidopsis mutants Genes enriched or selectively expressed in pollen and down-regulated in flower tissues (including pollen) of both ga1-3 and 35S:miR159a mutants are listed. "Pollen Selective" and "Pollen Enriched" data are from Pina et al. (2005). "Changed in Growth" indicates that gene expression is down-or up-regulated (at least 1.6-fold difference) during tube growth; data are from Wang et al. (2008). "Penta/Wild Type" shows a relative gene expression level in flower tissues (including pollen) of a penta mutant over wild-type Arabidopsis. The penta mutant carries five mutations in which gal1-3 and the members of the DELLA gene family are knocked out (Cao et al., 2006). "HSP90RNAi-A3/Wild Type" shows a relative gene expression level in flower tissues (including pollen) of a hsp90 knockdown line over wild-type Arabidopsis (Sangster et al., 2007). The data of "Penta/Wild Type" and "HSP90RNAi-A3/Wild Type" are from Genevestigator V3 (Hruz et al., 2008). "Reported" indicates that defective pollen tube growth in a knockout line was reported previously (Boavida et al., 2009;Szumlanski and Nielsen, 2009 . Numbers and percentages of vesicle trafficking genes that are present, enriched, or selectively expressed in Arabidopsis pollen. A total of 328 vesicle trafficking genes that are expressed at statistically high levels (detection call in the Affymetrix analysis is present) in pollen, sperm, seedling, silique, cotyledon, leaf, petiole, stem, root, or root hair zone were analyzed (Honys and Twell, 2004;Pina et al., 2005;Borges et al., 2008). Enriched indicates genes that express in pollen at least 1.2-fold higher than in other tissues (Honys and Twell, 2004;Pina et al., 2005;Borges et al., 2008). Selective indicates genes whose detection call is present in pollen but absent in the other tissues. GA regulated indicates genes whose expression is down-regulated in Arabidopsis mutants of the GA signaling pathway identified by Genevestigator V3 (Hruz et al., 2008). Forty-nine percent of vesicle trafficking genes are present and 51% are absent. Forty-two genes are enriched or selectively expressed in pollen. These data indicate that selected vesicle trafficking genes are used, and 14 of them are upregulated by the GA signaling pathway, in pollen. in seedlings of the HSP90RNAi-A3 line (Table II; Sangster et al., 2007). In animals, HSP90 is known to regulate Rab GTPase recycling in vesicle trafficking (Chen and Balch, 2006). Although the exact reason for up-regulation of these genes in seedlings is not yet clear, we speculate that when the HSP90 protein is expressed at a certain low level in seedlings, a positive feedback to the gene expression may occur to maintain normal vesicle trafficking. Among these 14 genes, AtRABA4d and AtAGD10 (Table II) knockout mutations already have been shown to impair pollen tube growth (Boavida et al., 2009;Szumlanski and Nielsen, 2009). Hence, we believe that this transcriptome analysis elucidated key vesicle trafficking genes that regulate pollen tube growth, even if not all vesicle trafficking genes were identified. This analysis also suggests that the network of molecules that regulates vesicle trafficking through the GA signaling pathway may be decoded by expanding the search to other gene families such as cytoskeletons, calcium channel proteins, and the GA signaling pathway in the future.

A Mathematical Model of Vesicle Trafficking in Pollen Tubes
To study the function of vesicle trafficking machinery at a systems level, we constructed a mathematical model that connects the functions of molecules expressed in pollen tubes, which includes gene products elucidated by the transcriptome analysis, using a set of ordinary differential equations (for details, see "Materials and Methods" and Supplemental Proof S1). Heinrich and Rapoport (2005) originally constructed a model to describe vesicle trafficking between the ER and Golgi apparatus. The model was used to computationally explain what kind of initial conditions satisfy equilibrium (generation of nonidentical compartments: the ER and Golgi apparatus) with limited molecular machinery (SNAREs and coat proteins). We adapted this model to describe vesicle trafficking among the Golgi apparatus, pollen tube tip, and recycling endosome (Samaj et al., 2006;Fig. 2). Because of the geometry of the vesicles and pollen tube, the number of vesicles required to elongate the plasma membrane and construct the cell wall is Figure 2. Schematic drawing of a mathematic model of vesicle trafficking in a pollen tube. A, Three-compartment model. We assume that the Golgi apparatus, tube tip, and recycling endosome (large circles and a semicircle) are major compartments that are involved in the transportation of newly synthesized cargos, pectin, and cargo phospholipids (gray) for pollen tube growth. Pectin and cargo phospholipids are de novo synthesized in the Golgi apparatus before tube elongation. These cargos are transported to the tube tip via vesicle (small circles) trafficking established in the pollen. At the tube tip, pectin and cargo phospholipids are incorporated into the plasma membrane and cell wall so that the tube elongates. The majority of SNARE X and U, Y and V, and Z and W localize in the Golgi apparatus, tube tip, and recycling endosome, respectively, when the vesicle trafficking is at equilibrium. GTPase A, B, C, D, E, and F initiate specific vesicle trafficking pathways, as indicated by thick arrows. thought to be different (Derksen et al., 1995;Bove et al., 2008). In this mathematical model, vesicles are considered as transporters that carry two cargos: (1) phospholipids to expand the plasma membrane, and (2) pectin to construct the cell wall (Supplemental Fig.  S1). When a vesicle fuses at the tube tip (exocytosis), phospholipids and pectin are thought to be integrated into the side of the tube where elongation occurs (Supplemental Fig. S2). In biological samples, phospholipids are also used for the vesicle trafficking system (Carr and Novick, 2000). In the model, phospholipids that are used for tube elongation (designated cargo phospholipids) and for the vesicle trafficking system are distinguished (Supplemental Fig. S1). Therefore, the number of vesicles required to expand the plasma membrane and construct the cell wall is equal but the number of cargo phospholipids and pectins that each vesicle carries is different (Supplemental Fig. S3, A and B). Furthermore, frequencies of exocytosis and endocytosis at the tube tip simply rely on the interaction between cognate v-SNAREs and t-SNAREs and the activities of small GTPase (Arf or Rab GTPase), similar to vesicle budding and fusion at the ER and Golgi apparatus (Supplemental Fig.  S3C).
Previously, Heinrich and Rapoport (2005) quantified the parameter values that mathematically satisfied the conditions to maintain nonidentical compartments. Hence, we applied these parameter values in our model (for details, see "Materials and Methods").

Equilibrium Status of Vesicle Trafficking in a Computed Pollen Tube
First, we mathematically proved that these three compartments (the Golgi apparatus, tube tip, and recycling endosome) can reach equilibrium (generation of nonidentical compartments) with the parameter values we used (Supplemental Proof S1). We subsequently applied the equations to the XPPAUT (for X-Windows Phase Plane plus Auto) program that solves differential equations based on our assumptions regarding pollen tube growth (Supplemental Program S1). The model was run for 600 time steps (t 0 # t # t 600 ) with a resolution of 0.2. The total relative area of compartments and vesicles is always 1 in the model (Supplemental Proof S1; Supplemental Program S1). The computational results indicated that vesicle trafficking reaches equilibrium at t 25 . At the equilibrium, the relative area of the Golgi apparatus, in which the vast majority of membranes and SNAREs of the vesicle trafficking system are accumulated at t 0 , is reduced from 0.96 (t 0 ) to 0.07 (t 25 ). At the same time, relative areas of the tube tip and recycling endosome increase from 0.005 (t 0 ) to 0.21 (t 25 ) and 0.14 (t 25 ), Figure 4. Computational results of changes in pollen tube length and amounts of pectin or cargo phospholipids as a function of time. A, Changes of pollen tube length. The maximum relative length of a pollen tube is 1. B, Changes of pectin amount in the Golgi apparatus, tube tip, and recycling endosome. A total relative amount of pectin is always 1. C, Changes of phospholipid amount in the Golgi apparatus, tube tip, and recycling endosome. The total relative amount of pectin is always 1. The results indicate that a pollen tube grows logarithmically due to logarithmic transport of pectin and nearly linear transport of cargo phospholipids to the tube tip. The total surface area of the three compartments and vesicles is always 1. The results indicate that the vesicle trafficking reaches equilibrium at t 25 when it is computed with 0.2 resolution. The results also indicate that the total area of vesicles in a pollen tube is larger than the compartments when the vesicle trafficking is at equilibrium. respectively (Fig. 3). The relative total area of vesicles increases from 0.03 at t 0 to 0.57 at t 25 (Fig. 3). This indicated that the majority of the surface area required for vesicle trafficking is in vesicles at equilibrium.

Computed Pollen Tube Growth Is Similar to Biological Observations
Next, we computed changes in the length of a pollen tube (Fig. 4A), amounts of pectin (Fig. 4B), and cargo phospholipids (Fig. 4C) in the Golgi apparatus, tube tip, and recycling endosome as a function of time. In the model, the maximum length of the pollen tube and total amount of pectin and cargo phospholipids are all set to 1. The computational growth curve of the pollen tube is logarithmic. The length reaches 0.924 at t 524 and hardly changes after this time point. This is due to a logarithmic transport of pectin and nearly linear transport of cargo phospholipids at the tube tip. Pectin in the recycling endosome increases until t 60 (the relative amount is 0.19) and gradually decreases thereafter. Meantime, cargo phospholipids in the recycling endosome rarely increase over time. Logarithmic growth of in vitro-grown pollen was reported previously (Boavida and McCormick, 2007), which supports the results of this model.
To change these relative numbers to absolute numbers, we used previously published data in which areas of vesicles and pollen tubes of Arabidopsis are measured by a transmission electron microscope (Table III; Ketelaar et al., 2008). According to the data, the diameter of a vesicle is 0.182 mm; hence, the surface area of a vesicle is calculated at 0.10 mm 2 . The average diameter of pollen tubes is 4.53 mm. We assumed that the shape of the tube tip is a hemisphere with a diameter of 4.53 mm. Thus, the calculated surface area of the tip tube is 32.23 mm 2 . Based on these measurements, we estimated that the total membrane area of the vesicle trafficking system among the three compartments is 153.48 mm 2 . At equilibrium, the areas of Golgi apparatus and recycling endosome are 10.74 and 21.49 mm 2 , respectively (Table IV). The surface distributions in the Golgi apparatus and recycling endosome in Arabidopsis pollen have not been obtained. To determine if the model quantitatively predicts these values, an additional series of experiments needs to be conducted. On the other hand, Ketelaar et al. (2008) estimated, based on electron micrographs, that 1,486 vesicles would be present near the tube tip (Table IV). The number of estimated vesicles is 877 at equilibrium. This suggests that the estimated number of the computed vesicles is underestimated. However, this is probably due to the difference in the measurements between the model and the biological sample. The biological sample contains multiple compartments such as the ER and the vacuole (Derksen et al., 2002;Ketelaar et al., 2008) that are not part of the systems model.
We also did comparisons with other previously published biological data in which the lengths of Arabidopsis pollen tubes are measured in vitro as a function of time. The maximum length of in vitrogrown pollen tubes varied from 250 to 1,000 mm, depending on conditions (Boavida and McCormick, 2007;Ketelaar et al., 2008;Szumlanski and Nielsen, 2009). Given these data, we assumed that the maximum length of pollen is 1,000 mm at 16 h when pollen is grown in vitro under optimum conditions (Boavida and McCormick, 2007). Because we computed time from t 0 to t 600, we assumed that a single time unit represents 1.6 min in the model. Based on these measurements and assumptions, we estimated that vesicle trafficking in a pollen tube reaches equilibrium within 40 min (1.6 min 3 25) after the initiation of the pollen tube tip (t 0 ). We also calculated the growth rate of the pollen tube as it grows from 10 to 100 mm. The rates were measured in in vitro-grown Arabidopsis pollen (Szumlanski and Nielsen, 2009). We found that the average growth rate of pollen tubes is 3.12 6 0.52 mm min 21 (from 2.8 to 5.9 mm min 21 ) in our model (Table IV). Ketelaar et al. (2008) and Szumlanski and Nielsen (2009) report that the growth rate of in vitro-grown Arabidopsis pollen is 4.07 6 0.36 and 4.5 6 1.0 mm min 21 , respectively. These numbers indicate a reasonable agreement between computational estimates and biological observations in pollen tube growth.  Ketelaar et al. (2008). **, Data are from Szumlanski and Nielsen (2009). ***, Estimated when the vesicle trafficking is at equilibrium. ****, Estimated when a computed tube length is from 10 to 100 mm. n/a, Not available.

Source Golgi Apparatus
Recycling Endosome
Tube Growth Rate mm 2 mm 2 mm min 21 Biological observation n/a n/a 1,486* 4.0* or 4.5** Computational estimate 10.74 21.49 877*** 2.8 to 5.9**** We then tested whether or not our model could correctly estimate vesicle trafficking in mutated pollen. The growth defect of in vitro-grown pollen tubes in the Arabidopsis raba4d knockout line was previously reported (Szumlanski and Nielsen, 2009). Because the closest yeast homologs of AtRabA4d, Ypt31/ 32, are known to initiate vesicle budding from the trans-Golgi network and direct the vesicles to the plasma membrane (Benli et al., 1996;Jedd et al., 1997), we used the raba4d knockout data to compare the function of GTPase B, which initiates vesicle transport from the Golgi apparatus to the tube tip in the mathematical model.
Although Arabidopsis expresses three AtRabA4d paralogous genes, AtRabA4a (71% amino acid identity), AtRabA4b (70% amino acid identity), and AtRabA4c (77% amino acid identity), which are most likely functionally redundant to AtRabA4d, our transcriptome analysis revealed that the three paralogous genes are absent in pollen (Supplemental Table S1). Hence, we hypothesized that AtRabA4d (RabA4d hereafter) is the sole source that functions as GTPase B in pollen. Accordingly, we changed the parameter values that represent the function of GTPase B to 0 (Supplemental Program S1) so that the computational results represent the raba4d null mutant. The computational results showed that the average growth rate of the pollen tubes is 0.84 6 0.34 mm min 21 when we measured the same time window as in wild-type pollen (Supplemental Program S1). This suggests that the growth rate would be reduced to 26% (0.84 of 3.12) of the wild-type pollen rate (Table V). This number reasonably agrees with the observation of the Arabidopsis raba4d knockout line, in which the pollen tube growth rate is reduced to about 18% of that of wild-type pollen (Szumlanski and Nielsen, 2009).
We also estimated that the length of in vitro-grown pollen tubes of wild-type Arabidopsis is 570 mm, while the length of raba4d knockout pollen is 270 mm (48% of wild-type Arabidopsis) at the same time point (Fig. 5; Table V). Szumlanski and Nielsen (2009) show that the average length of in vitro-grown pollen tubes of raba4d knockout Arabidopsis is about 45% of the wild-type Arabidopsis length after 24 h of incubation when the average length of wild-type Arabidopsis pollen is about 570 mm. Once again, these numbers indicate that the computational results reasonably agree with observations of the Arabidopsis raba4d knockout line. However, the computational results estimated that the length of raba4d knockout pollen would be 94% of that of wild-type pollen after 16 h of in vitro culture when the length of wild-type pollen reaches 925 mm (Fig. 5). The growth of raba4d knockout pollen is indistinguishable from that of wild-type pollen in vivo (Szumlanski and Nielsen, 2009). Therefore, we assume that the length of raba4d knockout pollen grown in vitro would not be correctly estimated by our current model. Alternatively, the length of in vitro-grown raba4d knockout pollen might become indistinguishable from that of wild-type pollen if pollen is grown under optimum conditions.

Computational Results Predict Localizations of Pectin and SYP125 in raba4d Mutant Pollen
We further predicted localizations of SYP125 in the raba4d mutant, which have not been investigated in biological samples yet. Our computational results predicted that localizations of SYP125, one of the vesicle machineries regulated by the GA signaling pathway in Arabidopsis pollen (Table II), in the raba4d knockout are little different from those in wild-type pollen (Fig. 6A). On the other hand, the amount of pectin transported in the raba4d knockout tube tip is reduced to 43% of wild-type pollen (Fig. 6B). This is due to higher recruitment rates of SNAREs to vesicles Figure 5. Computational estimations of pollen tube length in wild-type (WT) and raba4d knockout Arabidopsis. Changes of pollen tube length in micrometers were estimated as a function of time in hours. The estimations were made by, first, changing computational parameters so that a mathematical function of GTPase B (wb = 1, 1/kxb = 1/kub = 1/ kpb = 0.01, and 1/kyb = 1/kvb = 1/kzb = 1/kwb = 1) is halted (wb = 0 and 1/kxb = 1/kub = 1/kyb = 1/kvb = 1/kzb = 1/kwb = 1/kpb = 0) and, second, applying biological data of in vitro-grown wild-type Arabidopsis pollen tubes. The results indicate that when the length of wildtype tubes is 570 mm, the length of raba4d knockout pollen tubes is 270 mm (48% of the wild-type value). At 16 h, the length of raba4d knockout pollen tubes is 94% of the wild-type pollen. Table V. Data comparisons between biological observations and computational estimates in raba4d knockout pollen *, Data are from Szumlanski and Nielsen (2009). **, Estimated in the time window when the length of a computed pollen tube of wild-type Arabidopsis is from 10 to 100 mm. ***, Estimated when the length of a computed pollen tube of wild-type Arabidopsis is 570 mm.
In the model, SYP125 is transported into the tube tip soon after the initiation of the tube tip (t 0 ) through the recycling endosome (Supplemental Fig. S4). Decreased transportation of pectin in the raba4d knockout tube tip was already reported (Szumlanski and Nielsen, 2009). It should be possible to test the prediction that the transportation of SYP125 is little different from that of wild-type Arabidopsis by expressing GFP-SYP125 (Enami et al., 2009) in raba4d knockout Arabidopsis.
Computational Results Predict a Possible Phenotype of syp124 and syp125 Knockout Lines We also predicted possible phenotypes of syp125 mutant lines that are not reported yet. We assumed that SYP124 (85% amino acid identity), which is also expressed in pollen, is functionally redundant with SYP125, because these two SNAREs belong to the same Qa subfamily. SYP121 (62% amino acid identity), SYP122 (52% amino acid identity), and SYP123 (63% amino acid identity) also belong to the same Qa subfamily. However, our transcriptome analysis suggested that SYP121, SYP122, and SYP123 are depleted or absent in pollen (Supplemental Table S1). Accordingly, in this prediction, we changed the total amount of SNARE Y from 0.5 to 0.25 or 0 in the model. We assumed that 0.25 represents a syp125 or syp124 (single) knockout and 0 represents a syp124/syp125 (double) knockout mutant. SYP131 (44% amino acid identity), the other Qa subfamily gene, is also highly expressed in mature pollen. However, we did not consider SYP131 as the t-SNARE in the mathematical model, because SYP131 does not show polarized tip localization as SYP124 and SYP125 do (Enami et al., 2009) but rather homogeneous localization in the pollen plasma membrane (Enami et al., 2009). Our current model is limited to only vesicle trafficking to the tube tip, not to the entire plasma membrane.
The results predicted that the length of syp124 or syp125 mutant pollen would not be different (99%) from that of wild-type pollen (Fig. 7). However, syp124/syp125 (double) knockout mutant pollen would be shorter (66% of the wild-type length). This is due to nonlinear algebraic equations in the mathematic model (Supplemental Proof S1). Knockout analyses in these genes are required to test this prediction. Figure 7. Computational predictions of pollen tube growth in syp124 and syp125 knockout Arabidopsis. Changes of pollen tube length in micrometers were predicted as a function of time in hours. The predictions were made by, first, changing computational initial conditions of SNARE Y from YY1(0) = 0.49, YY2(0) = 0.005, and YY3(0) = 0.005 to YY1(0) = 0.245, YY2(0) = 0.0025, and YY3(0) = 0.0025 (representing single gene knockout) or YY1(0) = 0, YY2(0) = 0, and YY3 (0) = 0 (representing double gene knockout) and, second, applying biological data of in vitro-grown wild-type (WT) Arabidopsis pollen tubes. The computational results predict that the length of pollen tubes in single gene knockout syp124 or syp125 is 99% (922 or 925 mm) of the wild-type pollen when pollen is grown in vitro with the most optimized conditions. On the other hand, the length of pollen tubes in double genes knockout of syp124 and syp125 is 66% (614 mm) of the wild-type pollen (925 mm). Figure 6. Computational predictions of SYP125 localizations in wild-type (WT) and raba4d knockout Arabidopsis. A, Relative amounts of SYP125 that localized in the Golgi apparatus, tube tip, or recycling endosome in a wild-type pollen tube (when the length is 570 mm) or in a raba4d knockout pollen tube (when the length is 270 mm). The total relative amount of SYP125 in a pollen is 0.5. B, Relative amounts of pectin that localizes in the Golgi apparatus, tube tip, or recycling endosome in a wild-type pollen tube (when the length is 570 mm) or in a raba4d knockout pollen tube (when the length is 270 mm). The total relative amount of pectin in a pollen is 1. The computational results predict that the localization of SYP125 in raba4d knockout pollen is almost indistinguishable from that in wild-type pollen, although pectin localization is significantly different.

REMARKS
We constructed a basic systems model that connects molecular functions of vesicle trafficking genes during pollen tube growth of Arabidopsis. We identified 14 genes, 10 of which are small GTPases (and their effecters) or SNAREs, which are up-regulated by the GA signaling pathway during pollen development. The functions of these gene products are used to test or further develop the mathematical model constructed in this study. The computational results of this model, which is limited to connect the functions of small GTPases, SNAREs, and two cargo molecules (pectin and phospholipids), reasonably agree with observations of both wild-type and mutant Arabidopsis pollen so far.
We need to test knockout phenotypes of the 14 genes and identify subcellular localizations of their products to clarify the biological functions. The mathematical model reported here does not directly evaluate other vesicle machinery (such as exocyst, PI enzymes, and cytoskeletal proteins), other cargos (such as cellulose synthase), and other compartments (such as the ER, vacuoles, and non-tube-tip plasma membrane) that also should affect pollen tube growth. Moreover, the pollen-stigma interaction and the GA signaling pathway are not incorporated into the model. Therefore, even if our predictions reasonably agree with biological data, this model does not fully represent pollen tube growth or is oversimplified. Lastly, this model does not attempt to actually reflect any biophysical and biomechanical principles of pollen tube growth, such as the mechanics of the cell wall and turgor . While there are limitations to the model, it does connect molecular functions important for pollen tube growth and reasonably predicts the biology of the system. The model can be changed to reflect a variety of complications such as those described above.

Transcriptome Analyses
Supplemental materials of papers published by Pina et al. (2005), Honys and Twell (2004), and Borges et al. (2008) were used to identify the vesicle trafficking genes that are present in Arabidopsis (Arabidopsis thaliana) pollen. Excel files in the supplemental materials were imported to Microsoft Office Access to identify genes in pollen. The vesicle trafficking genes in the list (Supplemental Table S1) were selected based on annotations that were previously published in different papers (Elias et al., 2003;Sanderfoot and Raikhel, 2003;Vernoud et al., 2003;Lin et al., 2004;Yoshizawa et al., 2006;Paul and Frigerio, 2007;Anders and Jurgens, 2008). Genevestigator V3 (Hruz et al., 2008), with its default setting, was used to identify mutants in which a cluster of vesicle trafficking genes are down-regulated. The supplemental material of a paper published by (Wang et al., 2008) was used to analyze changes of expression levels of the selected genes during pollen tube growth.

Mathematical Modeling
The model was constructed based on a model constructed by Heinrich and Rapoport (2005). In our model, we made the following assumptions.
(2) Vesicle trafficking among these compartments requires v-SNAREs and t-SNAREs for vesicle fusion and small GTPase (Arf or Rab GTPase) for vesicle budding as the minimum machinery ( Fig. 2A).
(3) SNARE X and Y, U and V, and Z and W are native to compartments 1, 2, and 3, respectively. t-and v-SNAREs are not distinguishable (i.e. SNARE X and U can be v-or t-SNARE), but their interactions are limited within a pair for their native compartment (i.e. SNARE X interacts with U but not Y, V, Z, and W; Fig. 2A). The total amount of SNARE X, U, Y, V, Z, and W in the cell is 0.5, respectively.
(4) Small GTPase B and E, A and D, and C and F initiate vesicle budding from compartments 1, 2, and 3, respectively ( Fig. 2A). These small GTPases also direct the vesicles to a specific compartment via a protein interaction cascade (i.e. GTPase B initiates vesicle budding specifically from compartment 1 and directs the vesicles specifically to compartment 2).
(5) Budding rates are defined by the compartment area and the activity of a specific small GTPase. A direct pathway from the plasma membrane to the Golgi apparatus has not been confirmed in biological samples, although it may exist. Hence, the budding rate of vesicles that transfer from the tube tip to the Golgi apparatus is extremely low. Therefore, the equation has the form where wb1 denotes a budding rate of vesicles that transfer from the Golgi apparatus to the tube tip initiated by GTPase B, wb denotes a constant unique to GTPase B, and s1 denotes the area of the Golgi apparatus. Let wa be a constant unique to GTPase A that initiates vesicle budding from the tube tip and directs the vesicles to the Golgi apparatus. Then, wb = 1 and wa = 0.01.
(6) When SNAREs are recruited from a compartment to vesicles, SNAREs that are not native to the compartment are recruited more than those native to the compartment (i.e. when vesicles bud from the Golgi apparatus, SNARE Y, V, Z, and W are recruited to a vesicle at a 100-fold higher frequency than SNARE X and U. This difference was mathematically identified by Heinrich and Rapoport (2005), which allows the system to reach equilibrium. Consequently, the equation has the form where sxb denotes a saturation function for carrying SNARE X to vesicles that are initiated by GTPase B; x1, u1, y1, v1, z1, and w1 denote concentrations of SNARE X, U, V, Y, V, Z, and W at the Golgi apparatus, respectively; and kxb, kub, kyb, kvb, kzb, and kwb denote relative dissociation constants (reverse numbers of recruitment rates) of SNARE X, U, Y, V, Z, and W for GTPase B, respectively. kxb = kub = 100 and kyb = kvb = kzb = kwb = 1. (7) The rate of vesicle fusion depends on the concentration of interaction pairs of v-SNAREs and t-SNAREs in vesicles and a compartment. Therefore, the equation takes the form where fb2 denotes a fusion frequency of vesicles from the Golgi apparatus to the tube tip; K denotes a constant; xb1, ub1, yb1, vb1, zb1, and wb1 denote concentrations of SNARE X, U, Y, V, Z, and W in vesicles budded from the Golgi apparatus by GTPase B, respectively; and u2, x2, v2, y2, w2, and z2 denote concentrations of SNARE U, X, V, Y, W, and Z in the tube tip, respectively. K = 40.
(8) The vast majority of membranes and SNAREs that are required for the vesicle trafficking are accumulated in the Golgi apparatus at the initiation of a pollen tube (t = t 0 ; Fig. 2B). At a certain point in time (t = t n ), vesicle trafficking among the Golgi apparatus, tube tip, and recycling endosome reaches equilibrium. Consequently, initial conditions include where s1(0), s2(0), and s3(0) are relative compartment areas of the Golgi apparatus, tube tip, and recycling endosome at a pollen tube initiation (t = t 0 ), respectively. nb1(0), ne1(0), na2(0), nd2(0), nc3(0), and nf3(0) are total vesicle areas budded from the Golgi apparatus, tube tip, and recycling endosome by GTPases B, E, A, D, C, and F, respectively, at the initiation of a pollen tube (t = t 0 ), respectively. The total relative area of compartments and vesicles is 1.
(9) All pectin and cargo phospholipids that are required for pollen tube elongation are already synthesized in the Golgi apparatus at time 0 and delivered to the tube tip through vesicle trafficking as a function of time (Fig. 2B).
(10) The amount of pectin that vesicles can carry depends on a concentration of pectin in a compartment and a recruitment rate defined by a specific GTPase. Vesicles initiated by small GTPases at the tube tip do not carry pectin. Hence, the equation has the form spb1 ¼ p1 kpb 1 þ p1 kpb where spb1 denotes the saturation function for carrying pectin in vesicles initiated by GTPase B, p1 denotes the concentration of pectin in the Golgi apparatus, and kpb denotes the relative dissociation constant (the reverse number of the pectin recruitment rate) for GTPase B. Let 1/kpd be a pectin recruitment rate of GTPase D. Then, 1/kpb = 0.01 and 1/kpd = 0. The total amount of pectin in the cell is always 1.
(11) Phospholipids in the vesicle trafficking system and those to expand the plasma membrane are distinguished. Phospholipids to expand the plasma membrane are considered as cargos (designated cargo phospholipids; Supplemental Figs. S1 and S2). The amount of the cargo phospholipids that vesicles can carry depends on only the concentration in the compartment. Vesicles initiated at the tube tip do not carry the cargo phospholipids. Hence, the equation has the form where shb1 denotes a saturation function for carrying cargo phospholipids in vesicles initiated by GTPase B, buffer denotes a constant (0.03), and h1 denotes the concentration of cargo phospholipids in the Golgi apparatus. Let shd2 be a saturation function for carrying cargo phospholipids in vesicles initiated by GTPase D. Then, shd2 = 0. The total amount of cargo phospholipids in the cell is always 1.
(12) The tube length is defined by the amount of pectin and cargo phospholipids transported to the tube tip. Hence, the differential equation has the form where dl/dt denotes the rate of tube growth and mpb1, mpc3, mhb1, and mhc3 denote pectin and cargo phospholipids delivered to the tube tip from the Golgi apparatus and recycling endosome, respectively. The maximum length is 1.
Differential equations that describe our assumption were established and proved that the vesicle trafficking reaches equilibrium with the variables and constants we used (Supplemental Proof S1).

Computation
Codes in the ODE file made by Heinrich and Rapoport (2005) were modified based on our assumptions (Supplemental Program S1). The XPPAUT 5.98 program (http://www.math.pitt.edu/~bard/xpp/xpp.html/) was used to compute the differential equations. The results were exported to Microsoft Excel files, and the graphs were generated in the files.

Supplemental Data
The following materials are available in the online version of this article.
Supplemental Figure S1. Vesicle composition in the mathematical model. Figure S2. Exocytosis and endocytosis in the mathematical model.

Supplemental
Supplemental Figure S3. Predicted changes of the concentrations of pectin and cargo phospholipids in vesicles.
Supplemental Figure S4. Predicted changes of SYP125 and pectin amounts at the tube tip.
Supplemental Table S1. List of the vesicle trafficking genes.
Supplemental Proof S1. Mathematical proof of vesicle trafficking equilibrium.
Supplemental Program S1. Differential equations of vesicle trafficking. ACKNOWLEDGMENT