Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biotechnol Bioeng. Author manuscript; available in PMC 2013 February 1.
Published in final edited form as:
PMCID: PMC3240705

A population balance equation model of aggregation dynamics in Taxus suspension cell cultures


The nature of plant cells to grow as multicellular aggregates in suspension culture has profound effects on bioprocess performance. Recent advances in the measurement of plant cell aggregate size allow for routine process monitoring of this property. We have exploited this capability to develop a conceptual model to describe changes in the aggregate size distribution that are observed over the course of a Taxus cell suspension batch culture. We utilized the population balance equation framework to describe plant cell aggregates as a particulate system, accounting for the relevant phenomenological processes underlying aggregation, such as growth and breakage. We compared model predictions to experimental data to select appropriate kernel functions, and found that larger aggregates had a higher breakage rate, biomass was partitioned asymmetrically following a breakage event, and aggregates grew exponentially. Our model was then validated against several data sets with different initial aggregate size distributions and was able to quantitatively predict changes in total biomass and mean aggregate size, as well as actual size distributions. We proposed a breakage mechanism where a fraction of biomass was lost upon each breakage event, and demonstrated that even though smaller aggregates have been shown to produce more paclitaxel, an optimum breakage rate was predicted for maximum paclitaxel accumulation. We believe this is the first model to use a segregated, corpuscular approach to describe changes in the size distribution of plant cell aggregates, and represents an important first step in the design of rational strategies to control aggregation and optimize process performance.

Keywords: Population balance equations, plant cell culture, cell aggregation, bioprocess engineering, applied biological modeling, paclitaxel, Taxus


Plant cell culture is an important technology for the renewable supply of many natural products which are difficult to produce by traditional methods such as natural harvest or chemical synthesis. Though dedifferentiated plant cells can be grown in suspension culture and cultivated on a large scale using bioprocesses similar to other industrial cell cultures, the application of this technology is limited by low yields, slow growth rates, the tendency of plant cells to aggregate, and variability in all of these properties (Huang and McDonald 2009; Kolewe et al. 2008). In particular, plant cells grow as aggregates from two to thousands of cells (50 μm to well over 2000 μm in diameter). Cell aggregation has long been implicated in affecting the metabolism of individual cells (Kubek and Shuler 1978; Street et al. 1965) and many attempts have been made to understand and quantify the effect of aggregates on secondary metabolite production (Edahiro and Seki 2006; Hanagata et al. 1993; Hulst et al. 1989; Kessler et al. 1999; Kolewe et al. 2011; Madhusudhan and Ravishankar 1996; Pepin et al. 1999; Zhao et al. 2003).

Despite the importance of aggregation as a fundamental characteristic of plant cell cultures, a rigorous characterization of aggregation dynamics has remained elusive. Detailed information regarding aggregate size distributions has traditionally been difficult to obtain, as the prevalent measurement technique, mechanical sieving (Kessler et al. 1999; Mavituna and Park 1987; Zhao et al. 2003), provides low resolution data. As a result, most studies have tended towards qualitative and semi-quantitative observations. Early studies based on microscopy indicated that an increase in cell aggregate size was correlated with culture growth, and that a decrease in aggregate size was observed once growth stopped (Street et al. 1965; Wallner and Nevins 1973). Later studies confirmed the same behavior using size distributions obtained from sieves (Mavituna and Park 1987), and the idea of aggregates transitioning into different size classes over the course of batch subculture provided the basis for simple models to describe these changes (Hanagata et al. 1993; Kuboi and Yamada 1978). Various other techniques have been explored to measure aggregate size (McDonald et al. 2001; Pepin et al. 1999), but only recently has a reliable method been developed to enable the routine collection of aggregate size distributions (Kolewe et al. 2010). Aggregate size can now be monitored as any other process variable, enabling both a more comprehensive evaluation of process performance regarding metabolite production (Kolewe et al. 2011), as well as a more detailed mathematical description of aggregation dynamics, which is presented here.

A variety of models have been proposed to describe growth and product formation in plant cell cultures. Significant effort has been made to clarify the impact of nutrients on growth dynamics in plant suspension and hairy root cultures by developing models of increasing complexity (Cloutier et al. 2008; Curtis et al. 1991; Mairet et al. 2010; Vangulik et al. 1992). These models however are unsegregated and treat biomass as a culture-averaged parameter, and do not account for the heterogeneity of biomass induced by aggregation. The segregated models which have been developed typically divide biomass into general classifications (i.e., active and non-active) and use kinetic descriptions to account for growth dynamics (Bailey and Nicholson 1989; Schlatmann et al. 1999; Sirois et al. 2000; Vangulik et al. 1993). While providing some detail regarding the heterogeneity of biomass and indirectly accounting for aggregates of different sizes via the biomass classifications, these models do not explicate the phenomena underlying the transitions amongst classifications. A corpuscular approach has been used based on the cell cycle of individual plant cells (Degunst et al. 1990) but these results are difficult to compare to experimental data as plant cells do not grow as single cells. In this work we also take a corpuscular approach, but model plant cell aggregates as the basic cellular unit.

Population balance equations (PBEs) are used to model particulate systems in a variety of fields and applications ranging from solid-liquid systems such as crystallization to emulsions in food processing to cellular systems (Ramkrishna and Mahoney 2002; Ramkrishna 2000). PBEs are used to describe how a particle distribution changes over time based on kernel functions representing the relevant phenomena which affect the particles. In applied biology, PBEs have been most often used for systems with single cell particles such as yeast or bacteria (e.g., Henson 2003; Kromenaker and Srienc 1991; Mantzaris et al. 1999; Tsuchiya et al. 1966). They have also been used to describe multicellular systems such as filamentous fungi (Lin et al. 2008; Liu et al. 2005) and plant hairy root cultures (Han et al. 2004). There are no reports to our knowledge of PBEs used to describe aggregates in plant cell suspension culture, a problem which presents unique challenges including the lack of mechanistic detail available regarding plant cell aggregate growth and breakage processes.

In this work, we present a PBE model to describe and predict changes in the size distribution of plant cell aggregates that are observed experimentally. We select kernel functions and fit unknown kernel parameters based on a comparison of model predictions and experimental data. The model is then validated against several non-fitted experimental datasets with different starting distributions. Finally, the model is used to qualitatively evaluate the effect of culture shear environment by predicting paclitaxel accumulation from aggregate distributions corresponding to different breakage rates, demonstrating potential usefulness of the approach for developing process optimization strategies.

Materials and Methods


Taxus cuspidata cell line P991 suspension cell cultures were maintained in our laboratory as previously described (Kolewe et al. 2010), and plant cell aggregate size distributions were measured using a Coulter counter technique (Beckman Coulter, Inc., Brea, CA) (Kolewe et al. 2010). The measured aggregate volume was previously shown to correlate with biomass dry weight with ρ = 2.44 g DW / mLCoulter Counter volume (Kolewe et al. 2010), so all distribution data are presented here in terms of dry weight. For batch experiments, 200 mL working volume shake flasks were run in triplicate, and samples were taken every 2-3 days for aggregate size distribution measurements via the Coulter counter method and extracellular sugar measurements via a YSI 2700 Select Biochemistry Analyzer (YSI Life Sciences, Yellow Springs, OH). To initiate cultures with different initial aggregate size distributions from the same innoculum pool, aggregates were passed over a sterile 500 μm nylon mesh filter and the permeate was used to initiate cultures in triplicate with the same initial biomass concentration as unfiltered cultures, as described previously (Kolewe et al. 2011).

Model Formulation

The PBE results from a number balance on particles with respect to a given state variable, in this case aggregate volume, v. Assuming that plant cell aggregates can both grow and break apart, the governing PBE includes terms for the disappearance of aggregates due to the formation of larger aggregates by growth, the disappearance of aggregates due to the formation of smaller aggregates by breakage, and the formation of aggregates due to the breakage of larger aggregates, and is given by:


where n(v,t) is the continuous number density function and n(v,t)dv is the number of aggregates in size range v to v + dv at time t, g(v, S’) is the growth rate for aggregates of size v and effective intracellular concentration of total sugar S’, Γ(v) is the breakage frequency for aggregates of size v, and p(v,v’) is the partitioning function describing the distribution of daughter aggregates of size v resulting from the breakage of mother aggregates of size v’, assuming each breakage event results in two daughter aggregates. The partitioning function is subject to the following constraints:


ensuring the function is a proper probability density function and that the volume of daughter particles is equal to the volume of the mother particle from which the daughters were derived.

It is evident that biomass is lost due to cell death and biomass degradation once a plant cell suspension culture passes through its stationary phase (Bailey and Nicholson 1989; Sirois et al. 2000). As the PBE (1) is formulated to conserve volume, an additional sink term must be added to account for the biomass loss which is observed experimentally. The breakage of multicellular aggregates is not a natural process, but instead results from the shear stress that is part of the artificial in vitro suspension culture environment. Breakage then would not be expected to occur cleanly along cell-cell boundaries and instead results in cell lysis and the generation of cellular debris, both of which contribute to a total reduction in biomass. This phenomenon was described by adding an additional term to (1), which accounted for the disappearance of aggregates associated with partitioning via a parameter, b, representing the fraction of biomass which does not partition into daughter particles upon a breakage event. This parameter was assumed to be independent of aggregate size, so all breakage events result in the same percentage of biomass loss. The resulting PBE is similar to (1) but includes an additional parameter to account for the loss in biomass associated with the aggregate birth term:


where b = 0 corresponds to no biomass loss, in which case (3) reduces to (1), and b = 1 corresponds to complete biomass loss. To capture time course changes in a batch culture, (3) was coupled to a substrate balance which accounted for substrate depletion as a result of cell growth:


where S is the total extracellular sugar concentration and Y is a constant yield coefficient. The intracellular substrate concentration was calculated from the extracellular substrate concentration using a first-order filter:


This phenomenological equation represents a number of physiological processes including substrate uptake and conversion to usable sugars, where α is the rate constant for these lumped processes and describes how quickly cells respond to environmental changes (Zhu et al. 2000).

The determination of appropriate forms for the kernel functions Γ(v), p(v,v’), and g(v,S’) was emphasized since one major objective of this study was to gain insights into growth and breakage phenomenon underlying the formation of plant cell aggregate distributions. Consequently, different functional dependencies were evaluated. The breakage frequency, Γ(v), was considered to be either size independent or size dependent as follows:


where a is an adjustable parameter. Size independent breakage is often used for qualitative predictions (Vanni 2000) and was used to describe breakage of filamentous fungi (Lin et al. 2008). A power law dependence on particle diameter has also been used in a number of systems for better quantitative predictions (Vanni 2000) or where breakage is a result of shear (Barthelmes et al. 2003). For this case we limited the breakage frequency to be linearly dependent on diameter or linearly dependent on volume, as shown in (6). A generalized partitioning function, the extended Hill-Ng binary power-law product distribution (Diemer and Olson 2002), was used:


where B is the beta distribution, and q is a parameter greater than zero. This function covers a broad range of distributions from equal partitioning where both daughter aggregates are the same size (q → ∞), to a random or uniform distribution (q = 1), to an asymmetric U-shaped distribution in which there is a higher probability for one small and large aggregate (0 < q < 1), and to an erosion-type distribution, which describes the stripping of an infinitesimally small aggregate, (q → 0) (Fig. 1). The growth kernel, g(v,S’), was dependent on both aggregate size and effective substrate concentration, g(v,S’) = r(v)m(S’). Several forms of size dependence growth were considered:


where μmax is an adjustable parameter. The size independent expression, r(v) = μmax , has been used to describe the growth of single cell organisms (Nielson and Villadsen 1994; Zhu et al. 2000). The function where the rate of increase of the radius is constant, r(v) = μmaxv2/3, has been used to describe the growth of filamentous fungi (Nielson and Villadsen 1994), and corresponds to preferential growth of cells on the outer surface of an aggregate. The exponential growth expression, r(v) = μmaxv, which corresponds to the uniform growth of all cells in the aggregate, has also been used to describe filamentous fungi (Lin et al. 2008). Monod saturation kinetics was used for substrate dependence of the growth rate:


where Km is an adjustable parameter.

Figure 1
Extended Hill-Ng binary power-law product partitioning function, where the parameter q determines the shape of the distribution. The daughter aggregate distribution resulting from the breakage of an aggregate in three limiting cases are: q = 0.1, a sharp ...

Numerical Solution

The PBE was solved by discretizing the number distribution and employing a fixed-pivot technique (Kumar and Ramkrishna 1996). The spatial derivatives were approximated according to a method presented by David et al. (1991), and the resulting ODEs were solved using the MATLAB integration code ode45 (Mathworks, Natick, MA). This discretization method ensured that both total aggregate number and volume were conserved, which was deemed particularly important for this work as the total volume, or biomass, is a critical property. The accuracy of the numerical method was assessed using test cases of breakage only and growth only to confirm that total aggregate number and volume were in fact conserved.

To facilitate comparisons to model predictions, Coulter counter measurements were smoothed in MATLAB using a locally weighted scatterplot smoothing regression combined with a moving average. This histogram was then converted to a continuous number density, interpolated, and re-binned to produce a discretized distribution with variable spacing α v3 (uniform spacing in diameter) on a grid of 750 points, which was determined to be required for numerical solution convergence. As the minimum aggregate size which could be measured was larger than the smallest aggregates known to exist, the measured distribution was interpolated using cubic splines to n(vmin,t) = 0, where vmin is the minimum size of an aggregate, taken in this case to be a single cell with a diameter of 30 μm. To directly compare model predictions to experimental data, simulated distributions were re-binned according to the original histogram grid spacing, and a residual error, Ψ, was calculated based on the discretized volume distribution Nv= v N(v,t), as follows:


where Nv(vi, j) and N^v(vi,j) are the measured and predicted, respectively, values of the discretized aggregate distribution at aggregate size vi and experimental time point j, and the error is summed over N+1 distribution points and M experimental time points. This residual allowed for a direct comparison of goodness-of-fit between distributions which differed significantly in their magnitudes, and values of Ψ were used to judge the quality of the model predictions.


A representative dataset covering the range of aggregate sizes normally found in culture was used to structure and fit the model. Distributions shifted to larger aggregate sizes as the area under the aggregate curve increased, indicating that aggregate size increased concurrent with growth (Fig. 2a). Taking the first moment of the continuous number distribution (or equivalently the cumulative volume of the discretized volume distribution)


and converting to biomass, a growth curve providing an overview of the entire batch process was constructed (Fig. 2b). The growth curve could be divided into three general phases: an initial growth phase where substrate limitations are not present; a substrate-limited growth phase; and a stationary/death phase where growth is negligible and biomass is lost. To systematically identify appropriate kernel functions and parameters, these phases of the growth curve were simulated individually in an attempt to isolate the individual kernels, and then these pieces were integrated to simulate the entire growth curve.

Figure 2
Experimental data used to structure and parameterize the model. (a) Aggregate size distributions for Day 0 through Day 12, the period over which total biomass increased. (b) Biomass growth curve, illustrating phases of initial growth which was not substrate ...

Stationary phase: breakage frequency and partitioning kernels

After sugar was depleted and growth stopped, breakage was the only phenomena which needed to be considered. Therefore, this data was used to identify the breakage and partitioning kernel functions independent of growth. The aggregate size distributions shifted to smaller aggregate sizes, and the area under the aggregate curve decreased (Fig. 3). A simplified breakage-only model, where the growth term was removed from (3), was used to discriminate (6) and parameterize (7), starting with Day 12 as the initial distribution and Day 16 as the final distribution. The best fit for each form of the breakage kernel was determined, where a, the breakage rate constant, q, the parameter determining the daughter distribution function, and b, the parameter accounting for the loss of biomass upon breakage, were chosen by trial-and-error to give the minimum Ψ value. Two combinations of breakage and partitioning fit the experimental data well: size-independent breakage combined with the partitioning most resembling an equal distribution (q = 50, b = 0.3) and v1/3 dependent breakage combined with a steep U-shaped partitioning function (q = 0 .1, b = 0.15) (Fig. 3). Neither combination could be definitively selected based on this analysis, but breakage proportional to v was clearly worse than the other two combinations as indicated by both the residual error and a qualitative comparison. Therefore, this combination was not considered for subsequent analysis.

Figure 3
Effect of the breakage kernel function on predicted distributions for a breakage only model, where model input was the Day 12 experimental distribution. (a) Γ = a, (b) Γ = av1/3, and (c) Γ = av. In each case, the breakage rate ...

Initial growth phase: growth rate kernel

Keeping the breakage and partitioning kernel functions and parameters fixed for the two best combinations from the death phase analysis, the growth rate functional dependency was evaluated using data from the initial growth phase, Day 0 to Day 6. This analysis clearly showed that an exponential growth kernel provided the best fit to the experimental data (Fig. 4). While the kernel corresponding to growth of cells only on the aggregate surface yielded better results than size-independent growth for both breakage-partitioning combinations, this kernel did not show an increase in mean aggregate size as observed experimentally. As the exponential growthfunction was the only kernel which qualitatively fit the experimental data (Fig. 4c,f), the othergrowth kernels were not considered for further analysis. Even though the exponential growth kernel in combination with the size-dependent breakage kernel gave a slightly better quantitative fit than the combination with the size-independent breakage kernel, the improvement was small and not significant enough to specify the breakage kernel.

Figure 4
Effect of growth kernel size dependence on predicted distributions from the initial growth phase, Day 0 to Day 6, using the breakage combinations from the breakage-only model and fitting the specific growth rate, μmax, to minimize the residual ...

To further differentiate between the two breakage and partitioning kernel combinations we incorporated the intracellular substrate concentration, S’, using (5) to account for the lag phase before growth, ran the simulation for an additional two days to allow comparison with another experimental data point, and simultaneously refit parameters for each of the kernel functions. The size-dependent breakage combination resulted in approximately 20% improvement over the size-independent combination (Fig. 5). The size-independent breakage kernel significantly overestimated the contribution of large aggregates, and this error was more pronounced at later time points (Fig. 5b). While the size-dependent breakage combination underestimated the contribution of small aggregates, particularly at the later time points (Fig. 5a), this combination fit the experimental data reasonably well and was used for further simulations.

Figure 5
Simulation of growth from Day 0 through Day 8. A linear, first-order differential equation for the intracellular substrate concentration was added to account for the slight lag phase, and all kernel functions were simultaneously refit to minimize the ...

Entire batch: substrate-dependent growth

To simulate the time course of the entire batch culture, the additional parameters Y and Km in (4) and (9) were estimated from data over the entire 16 day experiment, and Table I lists the final model parameters and functions used for these simulations. Figure 6 shows comparisons of measured and simulated results for: aggregate distributions from Day 0 to Day 12, corresponding to points where biomass continued to increase (Fig. 6a); the biomass growth curve (Fig. 6b); the mean aggregate diameter (Fig. 6c); the total number of aggregates (Fig. 6d); and the amount of sugar in the extracellular media (Fig. 6e). While the simulated distributions began to diverge from the measured distributions in the latter stages of the growth curve, the calculated results for total biomass and mean aggregate size fit the experimental data extremely well (Fig. 6b,c). The comparison of the total number of aggregates indicates that the model underpredicted the total number as culture time progresses (Fig 6d). This is due to an underestimation of small aggregates, which make up a much larger percentage of the total aggregate number compared to their contribution to the total biomass. Extracellular sugar was described reasonably well by the model, although the initial rate of depletion was underestimated (Fig. 6e), which was likely due to over-simplification of the uptake and storage of sugars as a first order process.

Figure 6
Comparison of experimental data and model predictions for the entire batch process. (a) Aggregate distributions corresponding to the culture growth period, Day 0 to Day 12. (b) Total biomass growth curve. (c) Mean aggregate size. (d) Total number of aggregates. ...
Table 1
Nominal functions and parameter values used for Figures Figures66--99.

As the model was fit to the experimental data presented in Figure 6, it could be expected to agree reasonably well based on the number of adjustable parameters available. To test whether the model was sufficiently extensible to be used in a predictive capacity, we simulated several other data sets using the same kernel functions and model parameters as used in Figure 6. Aggregates from the same innoculum pool as used to generate Figure 6 were filtered to obtain a smaller initial distribution. The results indicate that our model can quantitatively predict the dynamic behavior of aggregate distributions that are significantly different from those used to fit the model (Fig. 7). In fact, the fit obtained in Figure 7 was slightly better than that in Figure 6 according to the Ψ measure, even though some qualitative behavior, such as a continuation of growth for several days past that observed experimentally (Fig. 7b), was slightly worse. To test whether the model was sufficiently robust to describe cultures which may differ in underlying properties other than the initial aggregate distribution, we used data from batch experiments that were carried out one year previously. Again, the model provided reasonable predictions, particularly during the growth phase for the most critical properties of total biomass and mean aggregate size (Fig. 8), indicating the broad extensibility of the model over the range of aggregate sizes encountered under normal processing conditions.

Figure 7
Comparison of experimental data and model predictions for the entire batch process. The model is compared to a non-fitted dataset with a significantly smaller initial aggregate distribution than that used for parameter estimation (Fig. 6). (a) Aggregate ...
Figure 8
Comparison of experimental data and model predictions for the entire batch process. The model is compared to a non-fitted dataset obtained one year earlier than that used for parameter estimation (Fig. 6). (a) Aggregate distributions corresponding to ...

Predicted effect of breakage rate on paclitaxel accumulation

To demonstrate the potential usefulness of the model and illustrate how it could be used in the development of bioreactor operating strategies, we performed simulations to quantitatively predict paclitaxel accumulation under different breakage rates, which could correspond to different shear environments induced by agitation. Though the relationship between paclitaxel accumulation and aggregate size is not simple (Kolewe et al. 2011), we modeled it using a simple sigmoidal curve, where smaller aggregates produce more paclitaxel, and paclitaxel production decreases to the limit zero as aggregate size increases (Fig. 9a). An additional ODE was used to couple paclitaxel accumulation to the PBE:


where T is the paclitaxel concentration and q(v) is the paclitaxel production rate. As smaller aggregates produce more paclitaxel (Kolewe et al. 2011), an obvious strategy to enhance paclitaxel accumulation would be to increase the number of small aggregates by increasing the breakage rate. However, an increase in breakage will also result in a greater loss of biomass according to our model. We performed a series of simulations starting from an experimentally determined aggregate distribution on Day 8, which is when cultures are normally elicited to stimulate paclitaxel production (Kolewe et al. 2011). The breakage rate constant, a, was varied to alter the effective breakage rate, and each simulation was run for two weeks post elicitation. The relative paclitaxel concentration is shown for a series of breakage rate constants in Fig. 9b. These results indicate that there is an optimum breakage rate to maximize paclitaxel accumulation, which occurred due to the balance of higher production in smaller aggregates with the loss of biomass at higher breakage rates.

Figure 9
Simulated paclitaxel accumulation under different shear conditions. (a) Sigmoidal representation of the dependence of paclitaxel production on aggregate size. (b) Total paclitaxel accumulation and mean aggregate size after two weeks in batch cultures ...


Kernel functions provide qualitative insight into aggregation phenomena

We have developed a population balance equation (PBE) model to predict the evolution of Taxus plant cell aggregate size distributions based on a phenomenological description of events which affect the aggregates. Rather than assume specific functional dependencies for the kernels representing the relevant phenomena from the outset, we evaluated a number of kernel functions that represented a range of realistic physiological possibilities, and selected the functions which best fit measured aggregate distributions from batch culture experiments. These phenomena can be difficult to measure independently, as breakage and partitioning, for instance, are inexorably linked, so this approach helps to provide insight into the qualitative nature of these processes. A size-independent breakage function has been previously used in cell-aggregating systems, though it was suggested that a size-dependent breakage function may be more appropriate (Lin et al. 2008). Though we did not account for shear effects explicitly, shear-dependent breakage has been modeled in other cellular systems using a size-dependence of v1/3 (Barthelmes et al. 2003). This description of breakage, while highly simplified from a mechanistic perspective, fit our data better than the size-independent function.

We found that the partitioning of aggregates after a breakage event was most accurately modeled with a slightly U-shaped distribution, where the most likely outcome was one small and one large aggregate. While more computationally demanding than a random distribution or a perfectly equisized distribution, the U-shaped function best fit the experimental data and was deemed to be the most physiologically plausible scenario. It was previously reported that cultures initiated with only large aggregates quickly developed a population of very small aggregates in tobacco culture (Kuboi and Yamada 1978), indicating that smaller pieces broke off the larger aggregates, rather than an alternative scenario of all aggregates breaking in half. We incorporated a loss of biomass associated with each breakage event to account for the biomass decrease observed in the stationary phase of batch culture. The nominal value for the biomass loss, corresponding to the parameter b, was found to be 20%. The effects of increased shear are most often associated with an impact on cellular metabolism (Gong et al. 2006), but these results suggest that an increased breakage rate (indirectly associated with increased shear in this study) would also lead to a direct loss of biomass, which may have a greater impact on culture performance.

By evaluating three different growth rate functions, we determined that the growth rate of individual aggregates appeared to be an exponential function of aggregate volume. This function implies that all cells within an aggregate are actively dividing, as opposed to the alternative growth functions considered such as size-independent growth or growth just on the aggregate surface. Exponential growth is probably not entirely mechanistically realistic, as diffusion limitations associated with larger aggregates (Curtis and Tuerk, 2006; Hulst et al. 1989) will cause metabolic changes for cells in the interior. However, the PBE framework presented here is general enough that a more detailed description of aggregate metabolism and growth could be incorporated into the model. While the model fits based on these empirically determined kernel functions were qualitatively correct, higher quality fits than those presented here have been achieved with PBE models of other particulate systems. We considered this study to be a preliminary effort that demonstrates the utility of PBE modeling of plant cell aggregation and suggest that incorporation of additional mechanistic detail into the kernel functions will likely improve model fits and predictions.

Utilization of the model as a predictive tool to guide operating strategies

The PBE model presented here can both reproduce experimental data from batch cell cultures, and predict the behavior of batch experiments for which the model was not fitted. Importantly, we presented model predictions of the aggregate distributions, as opposed to only lumped descriptions such as total biomass and substrate concentration that are presented for PBE models in other aggregating systems such as plant hairy roots (Han et al. 2004) or fungal hyphae (Liu et al. 2005). This distinction is critical as our purpose for creating a segregated model was to describe the distribution of biomass, which has been shown to impact critical process parameters such as paclitaxel accumulation (Kolewe et al 2011). Furthermore, as we have shown that paclitaxel accumulation is affected not only by the mean aggregate size but also by the shape of the distribution, the ability to predict the full aggregate distribution is an important feature of this model. In contrast to a PBE model based on more detailed mechanistic descriptions of aggregation and growth where only one set of experimental data was described (Lin et al. 2008), we developed a model based on a more simplified depiction of the relevant phenomena and demonstrated that it agreed well with non-fitted data sets. Thus, the model is appropriate for predicting process performance and running simulations to guide experimental investigations.

An important result obtained from simulation of different breakage rates was that there is an optimum breakage rate for the production of paclitaxel. While our analysis did not account for the potential negative impact of shear on cellular metabolism, this optimum was the result of a balance between the higher productivity of smaller aggregates and the loss of biomass associated with a higher breakage rate. These results point to the need for more sophisticated operating strategies aimed at managing aggregate size, as simply increasing aggregate breakage is not a viable approach. To this end, incorporating the effect of other process operating parameters, including shear and dissolved oxygen, as well as eventually including a more detailed structured description of nutrient utilization, would be highly desirable. Especially for plant cell cultures, where growth rates are very slow and culture performance is highly variable, simulations of culture conditions to optimize specific properties is highly preferable to large experimental arrays, which would be more suited to validation of the model predictions. This work represents an important step towards the development of strategies to control aggregate size, and the PBE model provides a basis to add additional details regarding both environmental conditions affected by reactor operating parameters as well as the state of individual cells within the aggregates.


This work was supported by grants from the National Science Foundation (CBET 0730779) and the National Institutes of Health (GM070852). M.E.K also acknowledges support from a National Research Service Award T32 GM08515 from the National Institutes of Health and the National Science Foundation-sponsored Institute for Cellular Engineering IGERT program (DGE-0654128).


  • Bailey CM, Nicholson H. A new structured model for plant cell culture. Biotechnol Bioeng. 1989;34(10):1331–1336. [PubMed]
  • Barthelmes G, Pratsinis SE, Buggisch H. Particle size distributions and viscosity of suspensions undergoing shear-induced coagulation and fragmentation. Chem Eng Sci. 2003;58(13):2893–2902.
  • Cloutier M, Bouchard-Marchand E, Perrier M, Jolicoeur M. A predictive nutritional model for plant cells and hairy roots. Biotechnol Bioeng. 2008;99(1):189–200. [PubMed]
  • Curtis WR, Hasegawa PM, Emery AH. Modeling linear and variable growth in phosphate limited suspension cultures of opium poppy. Biotechnol Bioeng. 1991;38(4):371–379. [PubMed]
  • Curtis WR, Tuerk AL. Oxygen transport in plant tissue culture systems. In: Gupta SD, Ibaraki Y, editors. Plant Tissue Culture Engineering. Springer; The Netherlands: 2006. pp. 173–186.
  • David R, Villermaux J, Marchal P, Klein JP. Crystallization and precipitation engineering .4. Kinetic-model of adipic acid crystallization. Chem Eng Sci. 1991;46(4):1129–1136.
  • Degunst MCM, Harkes PAA, Val J, Vanzwet WR, Libbenga KR. Modeling the growth of a batch culture of plant-cells - a corpuscular approach. Enzyme Microb Technol. 1990;12(1):61–71.
  • Diemer RB, Olson JH. A moment methodology for coagulation and breakage problems: Part 3 - generalized daughter distribution functions. Chem Eng Sci. 2002;57(19):4187–4198.
  • Edahiro J, Seki M. Phenylpropanoid metabolite supports cell aggregate formation in strawberry cell suspension culture. J Biosci Bioeng. 2006;102(1):8–13. [PubMed]
  • Gong YW, Li SY, Han RB, Yuan YJ. Age-related responses of suspension cultured Taxus cuspidata to hydrodynamic shear stress. Biochem Eng J. 2006;32(2):113–118.
  • Han BB, Linden JC, Gujarathi NP, Wickramasinghe SR. Population balance approach to modeling hairy root growth. Biotechnol Prog. 2004;20(3):872–879. [PubMed]
  • Hanagata N, Ito A, Uehara H, Asari F, Takeuchi T, Karube I. Behavior of cell aggregate of Carthamus tinctorius cultured cells and correlation with red pigment formation. J Biotechnol. 1993;30(3):259–269.
  • Henson MA. Dynamic modeling of microbial cell populations. Curr Opin Biotechnol. 2003;14(5):460–467. [PubMed]
  • Huang TK, McDonald KA. Bioreactor engineering for recombinant protein production in plant cell suspension cultures. Biochem Eng J. 2009;45(3):168–184.
  • Hulst AC, Meyer MMT, Breteler H, Tramper J. Effect of aggregate size in cell cultures of Tagetes patula on thiophene production and cell growth. Appl Microbiol Biotechnol. 1989;30(1):18–25.
  • Kessler M, ten Hoopen HJG, Furusaki S. The effect of the aggregate size on the production of ajmalicine and tryptamine in Catharanthus roseus suspension culture. Enzyme Microb Technol. 1999;24(5-6):308–315.
  • Kolewe ME, Gaurav V, Roberts SC. Pharmaceutically active natural product synthesis and supply via plant cell culture technology. Mol Pharmaceutics. 2008;5(2):243–256. [PubMed]
  • Kolewe ME, Henson MA, Roberts SC. Characterization of aggregate size in Taxus suspension cell culture. Plant Cell Rep. 2010;29(5):485–494. [PubMed]
  • Kolewe ME, Henson MA, Roberts SC. Increased accumulation of paclitaxel in Taxus suspension cultures based on the analysis of aggregate size as a process variable. Biotechnol Prog. 2011 In Press DOI: 10.1002/btpr.655. [PMC free article] [PubMed]
  • Kromenaker SJ, Srienc F. Cell cycle dependent protein accumulation by producer and nonproducer murine hybridoma cell lines - a population analysis. Biotechnol Bioeng. 1991;38(6):665–677. [PubMed]
  • Kubek DJ, Shuler ML. Generality of methods to obtain single cell plant suspension cultures. Can J Bot. 1978;56(20):2521–2527.
  • Kuboi T, Yamada Y. Changing cell aggregations and lignification in tobacco suspension cultures. Plant Cell Physiol. 1978;19(3):437–443.
  • Kumar S, Ramkrishna D. On the solution of population balance equations by discrettization. 1. A fixed pivot technique. Chem Eng Sci. 1996;51(8):1311–1332.
  • Lin PJ, Grimm LH, Wulkow M, Hempel DC, Krull R. Population balance modeling of the conidial aggregation of Aspergillus niger. Biotechnol Bioeng. 2008;99(2):341–350. [PubMed]
  • Liu G, Xing M, Han QG. A population-based morphologically structured model for hyphal growth and product formation in streptomycin fermentation. World J Microbiol Biotechnol. 2005;21(8-9):1329–1338.
  • Madhusudhan R, Ravishankar GA. Gradient of anthocyanin in cell aggregates of Daucus carota in suspension cultures. Biotechnol Lett. 1996;18(11):1253–1256.
  • Mairet F, Villon P, Boitel-Conti M, Shakourzadeh K. Modeling and Optimization of Hairy Root Growth in Fed-Batch Process. Biotechnol Prog. 2010;26(3):847–856. [PubMed]
  • Mantzaris NV, Liou JJ, Daoutidis P, Srienc F. Numerical solution of a mass structured cell population balance model in an environment of changing substrate concentration. J Biotechnol. 1999;71(1-3):157–174.
  • Mavituna F, Park JM. Size distribution of plant cell aggregates in batch culture. Chemical Engineering Journal and the Biochem Eng J. 1987;35(1):B9–B14.
  • McDonald KA, Jackman AP, Hurst S. Characterization of plant suspension cultures using the focused beam reflectance technique. Biotechnol Lett. 2001;23(4):317–324.
  • Nielson J, Villadsen J. Bioreaction Engineering Principles. Plenum Press; New York: 1994. p. 456.
  • Pepin MF, Smith MAL, Reid JF. Application of imaging tools to plant cell culture: Relationship between plant cell aggregation and flavonoid production. In Vitro Cell Dev Biol-Pl. 1999;35(4):290–295.
  • Ramkrishna D. Population balances: Theory and applications to particulate systems in engineering. Academic Press; San Diego: 2000. p. 355.
  • Ramkrishna D, Mahoney AW. Population balance modeling. Promise for the future. Chem Eng Sci. 2002;57(4):595–606.
  • Schlatmann JE, ten Hoopen HJG, Heijnen JJ. A simple structured model for maintenance, biomass formation, and ajmalicine production by nondividing Catharanthus roseus cells. Biotechnol Bioeng. 1999;66(3):147–157. [PubMed]
  • Sirois J, Perrier M, Archambault J. Development of a two-step segregated model for the optimization of plant cell growth. Control Eng Prac. 2000;8(7):813–820.
  • Street HE, Henshaw GG, Buiatti MC. Culture of isolated plant cells. Chem Ind. 1965;(1):27–33.
  • Tsuchiya HM, Frederickson AG, Aris R. Dynamics of microbial cell popultions. Adv Chem Eng. 1966;6:125–206.
  • Vangulik WM, Tenhoopen HJG, Heijnen JJ. Kinetics and stoichiometry of growth of plant cell cultures of Catharanthus roseus and Nicotiana tabacum in batch and continuous fermenters. Biotechnol Bioeng. 1992;40(8):863–874. [PubMed]
  • Vangulik WM, Tenhoopen HJG, Heijnen JJ. A structured model describing carbon and phosphate limited growth of Catharanthus roseus plant cell suspensions in batch and chemostat culture. Biotechnol Bioeng. 1993;41(8):771–780. [PubMed]
  • Vanni M. Approximate population balance equations for aggregation-breakage processes. J Colloid Interface Sci. 2000;221(2):143–160. [PubMed]
  • Wallner SJ, Nevins DJ. Formation and dissociation of cell aggregates in suspension cultures of pauls scarlet rose. Am J Bot. 1973;60(3):255–261.
  • Zhao D, Huang Y, Jin Z, Qu W, Lu D. Effect of aggregate size in cell cultures of Saussurea medusa on cell growth and jaceosidin production. Plant Cell Rep. 2003;21(11):1129–1133. [PubMed]
  • Zhu GY, Zamamiri A, Henson MA, Hjortso MA. Model predictive control of continuous yeast bioreactors using cell population balance models. Chem Eng Sci. 2000;55(24):6155–6167.