Effect of Spatial Segregation on Commensalistic ... - ACS Publications

Jun 2, 2011 - (11-13) Commensalism is an interaction in which one cell population (organism) is positively affected by the presence of the other, whil...
2 downloads 0 Views 1MB Size
ARTICLE pubs.acs.org/IECR

Effect of Spatial Segregation on Commensalistic Cultures—Series Reactors Satish J. Parulekar* Department of Chemical & Biological Engineering, Illinois Institute of Technology, Chicago, Illinois 60616, United States ABSTRACT: A comprehensive analysis of static and dynamic behavior of a mixed culture in two bioreactors in series is presented considering anaerobic digestion involving acidogens (X) and methanogens (Y) as the example bioprocess. A single continuous culture may operate at up to seven steady states, including up to four coexistence steady states, with only one coexistence steady state being locally stable. The bioreactors in the two-reactor system are identical in terms of composition of and dilution rate with respect to fresh feed. A two-reactor system may admit up to forty nine steady states, which are comprised of up to forty coexistence steady states, at least at very low interaction rates (R). When physically realizable, up to five coexistence steady states are locally stable. The static and dynamic analysis of the two-reactor system is facilitated by appropriate grouping of large number of steady states arising for very low R into nine clusters and grouping of the clusters into three cluster groups. Numerical illustrations reveal the rich steady state structure of the bioprocess for bioreactors in series. The two bioreactors can operate at up to five locally stable coexistence steady states over certain ranges of R. The two-reactor system is operationally more flexible and more robust vis-a-vis single reactor as concerns maintenance of mixed culture. Emergence of additional steady states at intermediate R reveals that series bioreactors are an example of a complex system.

1. INTRODUCTION Observed in their natural environment, most living species are found to exist in organized social settings reflecting various degrees of intermingling with other living species. The prospects of survival and coexistence of these species in these environments are strongly impacted by the nature of the interactions between them and the type of environment they populate. The ecosystem at large and any of its subsections are characterized by substantial spatial and temporal variations in key variables influencing the functioning of living species and responsible for biodiversity. Even in controlled settings in research laboratories and industrial complexes, these variations are common.1,2 Environmental partitioning in (bio)reactors arises when multiple reactors are used in series-parallel configuration or when multichamber reactors3 are used. Use of bioreactors (cell cultures) in series, considered here, is the norm for and an integral part of cultivation of living cells at bench, laboratory, pilot, and commercial scales.46 Living cell populations comprising of two or more different organisms and inhabiting a common environment are known to interact in a number of different ways.610 Mixed cultures involving wild-type and recombinant cell species are used extensively in commercial fermentations for production of a variety of food products, alcoholic beverages, and pharmaceuticals, biological waste treatment, and biological leaching, and are omnipresent in ecological systems.1113 Commensalism is an interaction in which one cell population (organism) is positively affected by the presence of the other, while the second population (organism) is not affected substantially by the presence of the first (see refs 6 and 1113 for citations on experimental and theoretical studies). The theoretical studies on commensalistic cultures have focused on operation of a continuous well-mixed bioreactor (continuous culture, CSTR) and have dealt with steady state characteristics11,12,1416 and dynamic behavior.11,12,14 The emergence phenomena to be r 2011 American Chemical Society

discussed later in this work has been demonstrated by periodic need-based addition of antibiotic to cultures containing recombinant and nonrecombinant cell species.17 Behavioral patterns of two coupled chemical reactors have been studied substantially, with the focus being on the study of chaotic dynamics, quasiperiodicity, and resonance phenomena resulting from the coupling of oscillators and synchronization of these.1820 This system can be viewed as an aggregation of two subsystems, the individual CSTRs, and displays the phenomenon of emergence featured by complexity theory. Emergence, in this sense, refers to any situation in which a system displays a level of functionality that is not possible for any of its subsystems when considered on their own. When the interaction between living species in a mixed culture is competition, the ability to support coexistence of these species is lacking from the homogeneous single reactor system,7 but emerges as a generic capability of the coupled reactor system with two-way exchange of cell culture between reactors.1,2124 In the present work, steady state and dynamic behavior of continuous mixed cultures of two cell species with commensalistic interaction and kinetic feedback in series bioreactors, identical as concerns composition of and dilution rate with respect to fresh feed, is analyzed and investigated. The specific example of commensalistic system considered in this study pertains to anaerobic digestion of complex insoluble organics. Anaerobic digestion is a biological process in which organic matter is transformed Special Issue: Nigam Issue Received: January 16, 2011 Accepted: June 2, 2011 Revised: May 29, 2011 Published: June 02, 2011 1525

dx.doi.org/10.1021/ie200103g | Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

Figure 2. Schematics of (a) a single well-mixed reactor and (b) two well-mixed reactors in series. J = A, B, C, X, Y.

Figure 1. Flowchart depicting steps involved in analysis of the tworeactor system. The interaction rate R is defined in eq 7.

by a community of microorganisms into biogas (primarily methane and carbon dioxide) in the absence of oxygen.2527 Because of its inherent advantages,28,29 this process is a promising method to solve some energy and ecological problems in agriculture and industry. In the past few years, a number of experimental studies aimed at improving the performance of anaerobic digestion30,31 and kinetic models for the process26,3236 have been reported. The process, involving a consortium of acid generating bacteria (host population) and methane generating bacteria (commensal population), requires hydolysis of complex insoluble organics by extracellular enzymes synthesized by the host population to produce simple soluble organics, which can be utilized by the host population. The flowchart in Figure 1 depicts the steps followed in analysis of series bioreactors. The problem formulation in section 2 is followed by steady state analysis for bioreactors in series in section 3, including (i) local stability analysis, (ii) classification, multiplicity, and stability of steady states for very weak one-way interaction between the two reactors, and (iii) variations in steady states with variation in one-way interaction. The essential background on single reactor operation14 is provided in the Appendix. Classification of steady states depends on whether or not host population and commensal population are retained in the two-reactor system. The two-reactor system can admit a large number of steady states. These steady states are organized into smaller number of clusters based on interaction between the host and commensal populations. Numerical illustrations on variations in steady state clusters and steady states and emergence of additional steady states with increasing interaction are provided in section 4. A qualitative discussion of dynamics of the two-reactor system follows in section 5. The one-way interaction between the two reactors implies that dynamics of the first reactor in the series is independent of that of the second reactor, while the dynamics of the latter is influenced by that of the former.

2. PROBLEM FORMULATION The anaerobic degradation of organic matter is a complicated biological process. The conversion of organic matter consists of

several independent, consecutive, and parallel reactions in which a close-knit community of bacteria cooperate to form a stable, self-regulating fermentation that transforms organic matter into a mixture of methane and carbon dioxide. These processes occur in six main stages, each carried out by one group of bacteria: hydrolysis of biopolymers into monomers (amino acids, sugars, and fatty acids), fermentation of amino acids and sugars, anaerobic oxidation of long chain fatty acids and alcohols, anaerobic oxidation of volatile fatty acids, conversion of acetate to methane, and conversion of hydrogen to methane.37 Many mechanistic models of varying complexity have been proposed for description of kinetics of anaerobic digestion.3740 While these models include multiple bacterial groups and many lumped chemical species and consider many reaction steps, a large number of kinetic parameters in these cannot be determined from the limited amount of experimental data available and must be assigned a priori. This places a serious limitation on predictive ability of these kinetic models. As an alternative, simpler kinetic models have been developed based on the complex mechanistic models to predict the dynamic behavior of anaerobic digesters. The six groups of bacteria are divided into two major groups: acid producing microorganisms or acidogens (X) and methane producing microorganisms or methanogens (Y).32,34,37,39 The simpler kinetic models also have a large number of parameters, which can be determined from the available experimental database.34 A concise three-stage representation of this complex process is14,32,34 E

X

Y

A f s B sf C f s P

ð1Þ

The fresh feed to each bioreactor comprises complex insoluble organics, A, which are hydrolyzed by extracellular enzymes, E, to simple soluble organics, B (Figure 2). The densities of various streams for a single reactor or a two-reactor system are considered to be the same, and the mass flow rates in and out of each reactor are considered to be equal. The extracellular enzymes are synthesized and secreted by the acidogens. The acidogens utilize B for their growth and generate volatile acids, C. The methanogens utilize C for their growth and generate, among other products, biogas, P, composed largely of methane and carbon dioxide. Background on the mechanism, kinetics, and modeling of anaerobic digestion is available in ref 14. The first step in eq 1 involves a feedback effect of X (E synthesized by X), adding to the complexity of the process. Such kinetic feedback also occurs in other cell cultures, two notable examples being conversion of cellulose and starch by suitable cell species to target metabolites. 1526

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

Utilization of these substrates requires synthesis of extracellular enzymes, cellulases in the case of cellulose and amylases in the case of starch, by living cells. The extracellular enzymes catalyze the hydrolysis of cellulose and starch to glucose, which can then be metabolized by living cells to allow for their growth and production of target metabolites. The conservation equations for two reactors in series, identical in terms of composition of and dilution rate with respect to fresh feed, are discussed next. The mass balances for reactor 1 are identical to those in eqs A1A5 with J = J01, J = A, B, C, X, Y, and D = D1, D1 = q1/V1. Since q1 = (q12 + q1e), (q2 + q12) = q2e, and (q1 + q2) = (q1e + q2e), from the total mass balances for the two bioreactors, the species conservation equations for the second bioreactor are 0

dA2 0 0 0 ¼ f12 , f12 ¼ D2 ðA0  A2 Þ  r12 + RðA1  A2 Þ dt

ð2Þ

0

dX2 0 0 0 0 ¼ f22 , f22 ¼  D2 X2 + ðμ12  k1 ÞX2 + RðX1  X2 Þ ð3Þ dt 0

dB2 0 0 0 0 ¼ f32 , f32 ¼  D2 B2 + r12  βμ12 X2 + RðB1  B2 Þ ð4Þ dt 0

dY2 0 0 0 0 ¼ f42 , f42 ¼  D2 Y2 + ðμ22  k2 ÞY2 + RðY1  Y2 Þ ð5Þ dt

Figure 3. Compartmentalization of the reaction system for (a) a single reactor and (b) two reactors in series.

All steady state solutions for a single reactor therefore are part of the portfolio of steady state solutions for the two-reactor system. Since the kinetics of generation/utilization of A, B, and X are not affected by concentrations of C and Y, both single reactor model and two-reactor model can be represented as twocompartment models (Figure 3). This compartmentalization simplifies considerably the analysis of the two-reactor system. Since q1e and q2e must be non-negative, R must satisfy the constraint V2 V1

D g ηR, η ¼

0

dC2 ¼ f52 , f52 dt 0

0

0

0

0

¼  D2 C2 + γμ12 X2  δμ22 Y2 + RðC1  C2 Þ

ð6Þ

with D2 ¼

q2 q12 , R ¼ V2 V2

ð7Þ

Serial culturing is a norm for cultivation of living cells with V2 being usually greater than or equal to V1 (V2 > V1 in scale-up operations). In such cases, R must be less than D.

3. STEADY STATE ANALYSIS OF REACTOR 2 The cell population balances for reactor 2 in eqs 3 and 5 reduce to

0

0

0

ð10Þ

0

0

ð11Þ

ðD + R + k1  μ12 ÞX2 ¼ RX1

The biogas production rate in reactor i, Qi, is Qi ¼ εμ2i Yi , i ¼ 1, 2

ð9Þ

ð8Þ

Because of the very restrictive online information, the control of the process is often reduced to regulation of the biogas production rate (Q). For this reason, Q is considered to be an important indicator of the performance of the mixed culture.27 The model and its special cases have been used in several prior studies and a large body of information on the ranges of the kinetic parameters is available.34 The robustness of the kinetic model, which captures the key regulatory effects in anaerobic digestion, has been tested in some of the prior studies. This was the motivation for the use of the model in this study. We consider here the special case D1 = D2 = D. The reactors need not have identical volumes and volumetric feed, effluent, and exchange rates, but are subject to the constraint q1/V1 = q2/V2. In subsequent discussion, R will be referred to as the interaction rate. The last term in each fk2, k = 15, accounts for one-way interaction between the two reactors. With the exception of the interaction term, the mass balances for reactor 2 in eqs 26 are identical in form to those for reactor 1 in eqs A1A5 and the former reduce to the latter as R f 0. For arbitrary R, if conditions in the two reactors are identical at a particular t, say t = 0, that is, if J01(0) = J02(0), J = A, B, C, X, and Y, then it follows from eqs 26 and A1A5 that the conditions in the two reactors will be identical in the transient operation for t > 0.

ðD + R + k2  μ22 ÞY2 ¼ RY1 X01

X02

One of the solutions of eq 10 is = = 0 and that of eq 11 is Y01 = Y02 = 0. Retention of X in the first reactor guarantees its retention in the second reactor with 0

0

ðD + R + k1  μ12 Þ > 0 for X2 > 0 if X1 > 0

ð12Þ

Similarly, retention of Y in the first reactor guarantees its retention in the second reactor with 0

0

ðD + R + k2  μ22 Þ > 0 for Y2 > 0 if Y1 > 0

ð13Þ

Retention of X in the second reactor when X01 = 0 requires that 0

0

ðD + R + k1  μ12 Þ ¼ 0 if X1 ¼ 0 and X2 6¼ 0 Similarly, retention of Y in the second reactor when requires that 0

0

ðD + R + k2  μ22 Þ ¼ 0 if Y1 ¼ 0 and Y2 6¼ 0

ð14Þ Y01

= 0 ð15Þ

The three types of steady states admissible for the two-reactor system are (a) total washout steady state (TW), X01 = X02 = 0, Y01 = Y02 = 0; (b) partial washout steady state (PW), Y01 = Y02 = 0 and X01 g 0, X02 > 0; and (c) coexistence steady state (CO), X01 g 0, X02 > 0, Y01 g 0, Y02 > 0. 1527

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

Table 2. Notation for Steady States and Stability Characteristics of These in a Two-Reactor System at Very Low Interaction (R f 0) Si, i = 17

single reactor

(1) TW, (2) PW1, (3) PW2, (4) CO1, (5) CO2, (6) CO3, (7) CO4 coupled reactors

Sij, i, j = 17. S = {Sij, i, j = 17} Sii (Sii  Si), i = 17, represent symmetric steady states.

Figure 4. Graphical representation of μ2 = (k2 + D). The curve indicates profile of μ2 and the horizontal line indicates (k2 + D). Cm = (K2KI)1/2, μ2m = μ20/(1 + (K2/KI)1/2)2.

total washout

S11 (stable for all R)

partial washout

S12, S13, S21, S22, S23, S31, S32, S33 (8 total)

coexistence

S1k, S2k, S3k, k = 47, S4j, S5j, S6j, S7j, j = 17 (40 total) steady states that are unstable for all R

Table 1. Possible Partial Washout and Coexistence Steady States steady state X (S2) PW1

C

Y

X1 C1 0

stability steady state X unstable

(S3) PW2

C

Y

X2 C2 0

coexistence

S2k, k = 47, S4j, S6j, S7j, j = 17 (25 total)

partial washout

S21, S22, S23 (3 total) steady states that are unstable for R f 0

stability stable/

coexistence

S14, S16, S17, S34, S36, S37, S52, S54, S12, S32 (2 total)

S56, S57 (10 total)

unstable (S4) CO1

X1 C1 Y1 unstable

(S5) CO2

X2 C1 Y2

stable/ unstable

partial washout

(S6) CO3

X1 C2 Y3 unstable

(S7) CO4

X2 C2 Y4

unstable

coexistence

S15, S35, S51, S53, S55 (5 total)

partial washout

S13, S31, S33 (3 total)

group 2 (G2)

steady states that are possibly stable for R f 0

group 3 (G3)

3.1. Number of Each Steady State Type for Very Low R. For very low R (R f 0), reactor 2 operates at a steady state that is in close vicinity of the steady state at which it will operate when not connected to reactor 1 (R = 0). If the two reactors operate at the same steady state type when uncoupled, they will operate at the same steady state type when coupled. The one-way interaction between the two reactors may lead in some cases to a change in the steady state type in reactor 2. Consider, for example, reactor 2 operating at total washout steady state and X being retained in reactor 1 at steady state. When the reactors are coupled, since B02 will be very low and μ12 = 0 (eq A7), the condition in eq 12 is satisfied and X02 > 0 (upgrade of steady state in reactor 2 from total washout to partial washout or coexistence). On a similar note, consider for example, reactor 2 operating at a total washout steady state and reactor 1 operating at a coexistence steady state. Retention of X in reactor 2 has been established above. When the reactors are coupled, since C02 will be very low and μ22 = 0 (eq A7), the condition in eq 13 is satisfied and Y02 > 0 (upgrade of steady state in reactor 2 from total washout to coexistence). If reactor 2 is operating at a partial washout steady state [C0 2 = C, C = C1, C2, Ci = γ(k1 + D)Xi/D, i = 1, 2], retention of Y in reactor 2 requires that the necessary condition in eq 13 be satisfied for reactor 2. This condition is of course violated for R f 0 when C1 (k2 + D) (Figure 4). Remark. Upgrading of steady state in reactor 2 from S2 to coexistence, when reactor 1 operates at a coexistence steady state (Sj, j = 47, Table 1), upon coupling is physically not realizable when C1 X1 > 0 (eq A9) in a single reactor operation, all nine clusters (Tables 2 and 3) are admissible. On the basis of the subsystem I characteristics for reactor 1, the nine clusters have been divided into three groups (Table 3): (A) C1, C7, and C8 (X01 = 0), (B) C2, C4, C9 (X01 = X1), and (C) C3, C5, C6 (X01 = X2) (X1 and X2 defined in eq A9). Operation of the two-reactor system at any steady state belonging to cluster group B is not possible due to instability of each steady state in it. The subsystem I state for reactor 1, A01, B01, and X01, is the same for each cluster and steady state in a cluster group and is invariant with respect to R. When reactor 1 operates at total washout steady state (S1) (cluster group A, A01 = A0, J01 = 0, J = B, C, X, Y), one can deduce that the mass balances for reactor 2 in eqs 26 reduce to those for a single reactor in eqs A1A5, with D being replaced by (D + R) and J = J02, J = A, B, C, X, Y, r1 = r12, μ1 = μ12, and μ2 = μ22 in the latter set (eqs 14 and 15). Clusters C7 and C8 are admissible for 0 < (D + R) < (μ1m  k1) and A0 g ψ(D + R) and merge with each other when A0 = ψ(D + R), with pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ψðD0 Þ ¼ B + ½βðk1 + D0 Þ + 4Rβðk1 + D0 ÞB, R K1 ðk1 + D0 Þ , D0 ¼ ðD + RÞ B¼ μ10  ðk1 + D0 Þ

λ7 ¼  D, λ8 ¼  ðD + RÞ, λ9 ¼ ðμ21  k2  DÞ, λ10 ¼ ðμ22  k2  D  RÞ Since these are real, cyclic states associated with a partial washout steady state will arise only if cyclic states are admissible for subsystem I, that is, only if AT admits a pair of imaginary eigenvalues. The limit points, characterized by λ9 = 0 or λ10 = 0, lead to emergence of coexistence steady states from particular partial washout steady states. For stability of a partial washout steady state, it is essential that max(λ9, λ10) < 0. If the two reactors have identical steady states, that is, if S0 1 = S02, S = A, B, C, X, and Y, then A1 = A2 and B1 = B2 in eqs 19 and

Admissibility of clusters C7 and C8 guarantees admissibility of partial washout steady states S12 and S13. The number of coexistence steady states in clusters C7 and C8 is obtained from Figure 5 with D being replaced by (D + R) as abscissa (region numbers indicated in parentheses): (1) S14, (2) S14, S15, (3) S15, (4) S15, S17, (5) S14, S15, S17, (6) S14, S16, (7) S14, S15, S16, S17, (8) S16, (9) S16, S17, and (10) S17. It follows from single reactor behavior14 that 1530

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research cluster C7 is unstable, cluster C8 is locally stable, steady states S12, S14, S16 and S17 are unstable, steady state S13 is unstable in regions 2 and 3 and locally stable in regions 4, 5, and 7, and steady state S15 is locally stable in regions 25 and 7 (portion of region 3 of interest). 3.3. Stability of Remaining Steady States at Very Low R. In a single reactor operation, steady states S2, S4, S6, and S7 are unstable.14 For the remaining steady states (S5j, j 6¼ 5, in regions 25 and 7 and S3j, j 6¼ 3, in regions 4, 5, and 7 of DA0 space), it should be evident from eq 20 that steady states S3j and S5j, j = 2, 4, 6, 7, if physically realizable, are unstable (Table 2). All steady states in clusters C6 and C7 (Table 3) are therefore unstable. Steady state S51 is locally stable in regions 25 and 7. In a single reactor operation, steady state S3 is locally stable in regions 4, 5, and 7 of the D-A0 space. Partial washout steady state S31 and coexistence steady states S35 and S53 are locally stable in regions 4, 5, and 7 (Table 4). While only one coexistence steady state, S5, is locally, asymptotically stable in regions 2, 4, 5, 7, and most of region 3 of the DA0 space for single reactor operation, the coupled reactor system permits operation at three stable coexistence steady states (S15, S51, and S55) in regions 2 and 3 and at five stable coexistence steady states (S15, S35, S51, S53, and S55) in regions 4, 5, and 7 and is therefore attractive and more flexible as concerns coexistence of X and Y. 3.4. Variations in Subsystems I and II for Bioreactor 2 with Variation in R. Extensive numerical effort has revealed that, with the exception of limit points, stability characteristics of the nine clusters and the physically realizable steady states belonging to these are unaltered upon variation in R. In regions 2, 4, 5, and 7 and portion of region 3 of interest (Figure 5), the locally, asymptotically stable clusters are C1 and C3 (symmetric) and C5 and C8 (asymmetric). Variation in R leads to variation only in steady states for reactor 2. At sufficiently high R, it is expected that J02 = J01, J = A, B, C, X, Y, and as a result, the two-reactor system will operate only at symmetric steady states. Satisfaction of the constraint in eq 9 at high R requires that V2/V1 be sufficiently small for operation of two-reactor system to be feasible. 3.4.1. Variations in Steady State Clusters. The asymmetric clusters (J01 6¼ J02, J = A, B, C) Cj, j = 49, will cease to exist beyond certain R as an asymmetric cluster merges with another asymmetric cluster, with at least one cluster in the pair being required to be unstable, at the limit points (|AT| = 0). Since each cluster contains a partial washout steady state, there is at least one admissible steady state in each cluster over the range of R where the cluster is admissible. For fixed D and A0, the number of clusters for the two-reactor system decreases from 9 (R f 0) to 3 (sufficiently high R). Owing to commonality of subsystem I characteristics for reactor 1 among clusters belonging to a cluster group, the merging cluster pairs are: C7C8, C4C9, and C5C6. 3.4.2. Variations in Asymmetric Steady States. For given D and A0, a particular asymmetric steady state Sij, if admissible, is a continuous function of R and will cease to exist beyond certain R as it will merge with another asymmetric steady state. For fixed D and A0, the number of steady states in the two-reactor system decreases from a maximum of n2 for R f 0 to n for sufficiently high R, with n being the number of steady states admissible for single reactor in a particular region in Figure 5 (Table 4). This reduction occurs through merges of certain asymmetric steady state pairs of one type (partial washout or coexistence), with at least one steady state in the pair being required to be unstable, at the limit points (eq 21).

ARTICLE

The merging steady states must correspond to identical conditions in reactor 1. Thus, asymmetric steady state Sij will merge with asymmetric steady state Sik, with i 6¼ j 6¼ k. Further, the merging steady states must belong to (i) a symmetric or asymmetric cluster or (ii) two merging asymmetric clusters. Where physically realizable and unless additional steady states emerge from these, the asymmetric steady state pairs subject to merger therefore are (1) Cluster C2 S24S26, S42S46, S62S64; (2) Cluster C3 S35S37, S53S57, S73S75; (3) Clusters C7 and C8 S12S13, S14S15 (D + R = Dc e D*), S16S17 (D + R = Dc e D*), S14S16 (D + R = D* e Dc), S15S17 (D + R = D* e Dc); (4) Clusters C4 and C9 S21S23, S25S27, S41S4i, S4jS4k, S61S6i, S6jS6k, i, j, k = 3,5,7, i 6¼ j 6¼ k; (5) Clusters C5 and C6 S31S32, S34S36, S51S5i, S5jS5k, S71S7i, S7jS7k, i, j, k = 2,4,6, i 6¼ j 6¼ k.

4. NUMERICAL RESULTS The steady states of a single reactor and coupled reactors for a specific set of operating parameters were identified using the continuation toolbox CL_MATCONT 1.141 for MATLAB 6.5. The software allows for tracing of steady state and periodic solutions, locating various bifurcation points, such as limit and Hopf points, and tracing continuation of bifurcation points. For the results to be presented next pertaining to variations in asymmetric clusters and steady states with variation in R, attributable to variations in J02, J = A, B, C, X, Y, the solid curves denote stable clusters or steady states, while the dashed curves denote unstable clusters or steady states. For this reason, only profiles of these are presented. Clusters C1, C2, and C3 being invariant with respect to R as concerns subsystem I for the tworeactor system, profiles of A02, B02, and X02 for these clusters are not shown. On a similar note, the symmetric steady states Sii, i = 17, being invariant with respect to R, profiles of A02, B02, C02, X02, and Y02 for these are not shown. The focus here is stable coexistence of X and Y in the two-reactor system. 4.1. Variations in Steady State Clusters. From the discussion at the end of section 3.4.2, it should be evident that clusters C7 and C8 are admissible for 0 < (D + R) e Dc and merge at R = (Dc  D) (limit point). As concerns clusters C4, C5, C6, and C9, there are two critical values of R, Rc1 and Rc2. The four clusters are admissible for R < Rc1. At R = Rc1, locally stable cluster C5 merges with unstable cluster C6. The other two clusters, C4 and C9, are unstable and admissible for 0 < R < Rc2 (Rc2 > Rc1) and merge at R = Rc2. None of the four asymmetric clusters are admissible for R > Rc2. The first illustration pertaining to variations in these four asymmetric clusters (Table 3) is for D = 0.29 day1 and A0 = 746 mg/L (Figure 6). The two critical R are Rc1 = 2.6559343  103 day1 and Rc2 = 2.1842324  102 day1. The unstable cluster C9 collides into unstable symmetric cluster C2 at R = 2.1749575  102 day1 (branch point B). The second illustration pertains to a nearby and higher D, D = 0.3 day1, and the same A0 (Figure 7). The two critical R are Rc1 = 7.40106  103 day1 and Rc2 = 8.3199036  103 day1. A comparison of the results in Figures 6 and 7 reveals that the limit points for cluster pairs C4C9 and C5C6 approach each other as D is increased. The third illustration pertains to a higher A0, A0 = 800 mg/L, and the same D (= 0.3 day1) (Figure 8). The two critical R are Rc1 = 3.203033  103 day1 and Rc2 = 1.6779186  102 day1. The unstable cluster C9 collides into unstable symmetric cluster C2 at R = 1.601715  102 day1 (branch point B). A comparison of 1531

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

Figure 6. Portraits of (a) R and A02, (b) R and B02, and (c) R and X02 for D = 0.29 day1 and A0 = 746 mg/L. The steady state composition for single reactor operation is reported in column 1 of Table 5. Clusters (i) C5 and C6 and (ii) C4 and C9 merge at the limit points L and cluster C9 collides into cluster C2 at the branch point B.

the results in Figures 7 and 8 reveals that the limit points for cluster pairs C4C9 and C5C6 approach each other as A0 is decreased. An insight into convergence of the two limit points with increase in D for fixed A0 or decrease in A0 for fixed D is obtained from the limit point continuation curves displayed in Figure 9. A cusp point occurs at A0 = ψ(D), where X1 = X2 (eq A9) and corresponds to collision of clusters (i) C4 and C5 and (ii) C6 and C9. The limit points for cluster pairs C4C9 and C5C6 diverge on a limit point continuation curve as one moves away from the cusp point. The coexistence steady state S51, which belongs to cluster C5, is locally stable for 0 < R < Rc1 and undergoes static bifurcation at R = Rc1 (represented by the lower portion of the limit point continuation curves in Figures 9a and 9b) to a steady state S5i, i = 2, 4, or 6, belonging to cluster C6. 4.2. Variations in Asymmetric Steady States. Since locally stable coexistence steady states belong to clusters C3, C5, and C8, illustrations are provided only for steady states belonging to these clusters. 4.2.1. Steady States in Cluster C3. 4.2.1.1. Regions 4, 5, and 7. Cluster C3 is comprised of nine steady states, at least for R f 0 (Tables 3 and 4), eight of these allowing for coexistence of X and Y. Steady states S73, S75 and S77 are unstable. Of the nine steady states, three (S33, S55, S77) are symmetric. There are three critical values of R, Rc3, Rc4, and Rc5, with Rc5 g Rc3 and Rc5 g Rc4. These represent limit points for coexistence steady state pairs (i) S35S37, (ii) S53S57, and (iii) S73S75, respectively (Figures 10 and 11). Steady states in a particular pair exist

ARTICLE

Figure 7. Portraits of (a) R and A02, (b) R and B02, and (c) R and X02 for D = 0.3 day1 and A0 = 746 mg/L. The steady state composition for single reactor operation is reported in column 2 of Table 5. Clusters (i) C5 and C6 and (ii) C4 and C9 merge at the limit points L.

over the same range of R. No asymmetric steady states in cluster C3 are admissible for R > Rc5. Let Rl = min(Rc3,Rc4) and Ru = max(Rc3,Rc4). When Rc3 6¼ Rc4, the number of steady states admissible in cluster C3 is (i) 9 for R < Rl, (ii) 8 for R = Rl, (iii) 7 for Rl < R < Ru, (iv) 6 for R = Ru, (v) 5 for Ru < R < Rc5, (vi) 4 for R = Rc5, and (viii) 3 for R > Rc5. When Rc3 6¼ Rc4, the number of locally, asymptotically stable steady states in cluster C3 is (i) 4 for R < Rl, (ii) 3 for Rl e R < Ru, and (iii) 2 for R g Ru. When Rc3 = Rc4, the number of steady states admissible in cluster C3 is (i) 9 for R < Rc3, (ii) 7 for R = Rc3, (iii) 5 for Rc3 < R < Rc5, (iv) 4 for R = Rc5, and (v) 3 for R > Rc5. When Rc3 = Rc4, the number of locally, asymptotically stable steady states in cluster C3 is (i) 4 for R < Rc3 and (ii) 2 for R g Rc3. One of these stable steady states, S33, is unwanted. The two illustrations provided next pertaining to variation in steady-state subsystem II composition for reactor 2 are for A0 = 1000 mg/L and D = 0.2 and 0.29 day1 (Figures 10 and 11). For D = 0.2 day1, the critical R values are Rc3 = 0.103799 day1, Rc4 = 2.567  103 day1, and Rc5 = 0.1489886 day1, and for D = 0.29 day1, these are Rc3 = 0.013799 day1, Rc4 = 0.0271433 day1, and Rc5 = 0.082634 day1. For fixed A0, Rc3 and Rc4 are expected to vary with D and there will be crossover of the profiles of Rc3 and Rc4 at a D in the interval 0.2 < D < 0.29. An insight into this is obtained from the limit point continuation curves for the six asymmetric steady states displayed for A0 = 1000 mg/L (A0 > A*) 0 in Figure 12a. The profiles of Rc3 and Rc4 exhibit a crossover. There are two cusp points for the limit point continuation curve for the steady state pair S73S75. The cusp point at higher D, D = D* = 0.304 day1, corresponds to merger of limit points for 1532

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

Figure 8. Portraits of (a) R and A02, (b) R and B02, and (c) R and X02 for D = 0.3 day1 and A0 = 800 mg/L. The steady state composition for single reactor operation is reported in column 3 of Table 5. Clusters (i) C5 and C6 and (ii) C4 and C9 merge at the limit points L and cluster C9 collides into cluster C2 at the branch point B.

steady state pairs S53S57 and S73S75, in agreement with merger of single reactor coexistence steady states S5 and S7. Since steady states S35 and S37 are admissible for 0 < (D + R) < D*, the limit point continuation curve for the steady state pair S35S37, (D + R) = D*, is a straight line (Figure 12a). The cusp point at the lower D, D = Dc1, Dc1 = 0.1642152 day1, corresponds to vanishing of steady state S7 (Y4 = 0) and merger of limit points for steady state pairs S35S37 and S73S75. Coexistence steady state S7 is not admissible and partial washout steady state S3 is unstable for D e Dc1. Based on the remark in section 4.1, steady state S53 is not admissible for D e Dc1 and as a result Rc4 f 0 as D f Dc1 (Figure 12a). This cusp point corresponds to transition from region 3 to region 4 with increasing D (Figure 5a). The continuation curve for limit points of steady state S35, (D + R) = D*, also applies for the interval 0 < D e Dc1, with steady states S3 and S35 being unstable. The limit point continuation curves for the six asymmetric steady states for a lower A0, A0 = 750 mg/L (A0 < A0*), in Figure 12b display patterns very different from those observed in Figure 12a. The profiles of Rc3 and Rc4 do not exhibit a crossover. Since steady states S35 and S37 are admissible for 0 < (D + R) < D*, the limit point continuation curve for steady state pair S35S37, (D + R) = D*, is a straight line (Figure 12b). The two cusp points at D = Dc1, Dc2, Dc1 = 0.207458 day1, Dc2 = 0.3001972 day1, correspond to vanishing of steady state S7 (Y4 = 0). Steady state S7 is not admissible and steady state S3 is unstable for D e Dc1 and D g Dc2. Based on the remark in section 4.1, steady state S53 is not admissible for D e Dc1 and

ARTICLE

Figure 9. Variation in R with variation in A0 (part a, D = 0.3 day1) or D (parts b and c, A0 = 746 mg/L) on the continuation curve for limit points of clusters (i) C4 and C9 (upper portions in parts a and b) and (ii) C5 and C6 (lower portions in parts a and b). A close-up of the lower portion of the limit point continuation curve in b is provided in c. C denotes the cusp point. The upper portions of the limit point continuation curves are the profiles of Rc2 and the lower portions of the same the profiles of Rc1.

D g Dc2 and as a result Rc4 f 0 as D approaches Dc1 and Dc2 (Figure 12b). The cusp point at lower D corresponds to transition from region 3 to region 4 with increasing D (Figure 5a) and the cusp point at upper D corresponds to transition from region 5 to region 2 with increasing D (Figure 5b). The cusp point associated with transition from region 5 to region 2 with increasing D is admissible only for A01 < A0 < A02 [A01 = ψ(D1), A02 = ψ(D2), D1 = 0.2783 day1, D2 = 0.3025 day1, Figure 5, parts b and c]. The continuation curve for limit points of steady state S35, line (D + R) = D*, also applies for the intervals 0 < D e Dc1 and Dc2 e D e Dc, with steady states S3 and S35 being unstable. For A02 e A0 < A*, 0 the limit point continuation curve for steady state pair S73S75 has only one cusp point, that associated with transition from region 3 to region 4 (Figure 5a). The limit point continuation curves for the six asymmetric steady states for D = 0.2, 0.25, and 0.29 day1 are displayed in Figure 13. The profiles of Rc3 and Rc4 cross each other at a critical A0, A0 = A0 (Rc3 > Rc4 for A0 < A0 and Rc3 < Rc4 for A0 > A0, Figure 13c). It is obvious from the profiles in Figure 13 that A0 decreases with increasing D. Since steady states S35 and S37 are admissible for 0 < (D + R) < D*, the limit point continuation curve for steady state pair S35S37, (D + R) = D*, is a horizontal line (Figure 13). The single cusp point at A0 = A0c, (i) A0c = 781.55 mg/L for D = 0.2 day1, (ii) A0c = 648.35 mg/L for D = 0.25 day1, and (iii) A0c = 696.35 mg/L for D = 0.29 day1, corresponds to vanishing of steady state S7 (Y4 = 0). Steady state S7 is not admissible and steady state S3 is unstable for A0 e A0c. 1533

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

Figure 10. Portraits of (a) R and C02 and (b) R and Y02 for steady states belonging to cluster C3 with D = 0.29 day1 and A0 = 1000 mg/L. Legend ij represents steady state Sij. The steady state composition for single reactor operation is reported in column 5 of Table 5. Steady states (i) S35 and S37, (ii) S53 and S57, and (iii) S73 and S75 merge at appropriate limit points L and steady state S73 collides into symmetric steady state S77 at the branch point B.

Based on the remark in section 4.1, steady state S53 is not admissible for A0 e A0c and as a result Rc4 f 0 as A0 f A0c (Figure 13). The cusp point corresponds to transition from region 4 to region 3 for D < D1, from region 4 to region 2 for D = D1 and from region 5 to region 2 for D1 < D < D2 as A0 is decreased (D1 = 0.2783 day1, D2 = 0.3025 day1, Figure 5, parts a and b). The continuation curve for limit points of steady state S35, the horizontal line R = (D*  D) also applies for the interval ψ(D) < A0 e A0c, with steady states S3 and S35 being unstable (Figure 13b). Regions 2 and 3 do not extend to the interval D2 e D < D*; the limit point continuation curves for variable A0 do not have cusp points in this interval as a result (Figure 5c). 4.2.1.2. Regions 2 and 3. Since only one coexistence steady state, S5, is admissible and steady state S3 is unstable in a single reactor operation in regions 2 and 3 of DA0 space, the stable cluster C3 in these regions is comprised of three steady states, S33, S35, and S55, with only steady state S55 being locally stable (Tables 3 and 4). The limit points for steady state S35 lie on the line (D + R) = D* (Figures 12 and 13). Another illustration for this is provided in Figure 14 for A0 = 600 mg/L, for which Du = 0.2805 day1 (eq A10). A cusp point is detected at D = Du, where A0 = ψ(D). Retention of X and therefore Y is not possible for D > Du. For 0 < D < Du, both single reactor and two-reactor system operate in regions 2 and 3 (Figure 5a). The unstable coexistence steady state S35 (Y02 > 0), which is admissible for R < (D*  D), cannot merge with symmetric steady state S55 as it is stable and cannot merge with partial washout steady state S33. At the limit point, steady state S35 must branch to an emerging coexistence steady state E, which is not admissible for R f 0. The profiles of C0 2 and Y0 2 in Figures 15a and 15b, for D = 0.25 day1 and A0 = 600 mg/L (region 3, Figure 5a) indeed reveal this to be the case. Steady state S35 is admissible up to R = 0.054 day1 and undergoes bifurcation to an unstable emerging coexistence steady state E. Steady state E collides with unstable steady state S33 at the branch point B (Figures 15a and 15b). The evolution of the emerging steady state E can be explained with profiles of Y02

ARTICLE

Figure 11. Portraits of (a) R and C02 and (b) R and Y02 for steady states belonging to cluster C3 with D = 0.2 day1 and A0 = 1000 mg/L. Legend ij represents steady state Sij. The steady state composition for single reactor operation is reported in column 4 of Table 5. Steady states (i) S35 and S37, (ii) S53 and S57, and (iii) S73 and S75 merge at appropriate limit points L.

Figure 12. Variation in R with variation in D on the continuation curves for limit points of asymmetric coexistence steady states belonging to cluster C3 with (a) A0 = 1000 mg/L and (b) A0 = 750 mg/L. C denotes a cusp point. The profiles of Rc3, Rc4, and Rc5 represent limit points for steady state pairs (i) S35S37, (ii) S53S57, and (iii) S73S75, respectively.

in Figure 15c for D = 0.25 day1 and A0 = 600, 648.34, and 750 mg/L. The reactor operation is in region 4 for A0 = 750 mg/L and on the interface between regions 3 and 4 for A0 = 648.34 mg/L. For both A0 values, steady state S35 is locally stable and bifurcates to unstable steady state S37 at the limit point. The profile for A0 = 648.34 mg/L is a nice illustration of upgrading of reactor 2 steady state from partial washout, S3, at R = 0 to coexistence, S7, for R > 0. The one-way interaction between the two reactors leads to upgrade in steady state from S33 to S37. For A0 < 648.34 mg/L, this upgrading is postponed until a certain R, when coexistence 1534

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

Figure 13. Variation in R with variation in A0 on the continuation curves for limit points of asymmetric coexistence steady states belonging to cluster C3 with (a) D = 0.2 day1, (b) D = 0.25 day1, and (c) D = 0.29 day1, C denotes a cusp point. The profiles of Rc3, Rc4, and Rc5 represent limit points for steady state pairs (i) S35S37, (ii) S53S57, and (iii) S73S75, respectively. The dashed line in (b) corresponds to unstable steady state S35, with A0 = ψ(D) at the lower cusp point.

Figure 14. Limit point continuation curve for coexistence steady state S35 with A0 = 600 mg/L. C denotes the cusp point corresponding to A0 = ψ(D).

steady state E emerges from steady state S33 (Figure 15). The emergence of steady state E, which is admissible in regions 2 and 3 and is unstable, leaves the number of locally stable steady states unaltered and dynamics of the two-reactor system unaffected. The emergence of coexistence steady state E is due to transition from simple (stand-alone operation of each reactor) to complex (interaction between two reactors). 4.2.2. Steady States in Clusters C5 and C8. Coexistence steady states belonging to clusters C7 and C8 are admissible for 0 < (D + R) < Du (eq A10) or 0 < R e Rc6, with Rc6 = (Du  D). The coexistence steady state pairs subject to merger at R = Rc6 are (1) S14S15 and S16S17 when Dc < D*, and (2) S14S16 and

ARTICLE

Figure 15. Portraits of (a) R and C02 and (b, c) R and Y02 for steady states S33, S35, and S37 and emerging coexistence steady state E, as appropriate, for D = 0.25 day1 and A0 = 600 mg/L (parts a, b and c) and 648.34 and 750 mg/L (part c). Legend ij represents steady state Sij. The steady state composition for single reactor operation at A0 = 600 mg/L is reported in column 6 of Table 5. Coexistence steady state S37 or E collides into partial washout steady state S33 at the branch point B. In c, legends A, C, and D represent profiles for A0 = 750, 648.34, and 600 mg/L, respectively.

S15S17 when Dc > D*. Steady states S31, S51, and S71 belonging to cluster C5, if physically realizable at R = 0, are admissible where this cluster is admissible (0 < R < Rc1). Each steady state in cluster C5 merges with an appropriate steady state of the same type in cluster C6, the merging steady state pairs being S31S32, S51S5i, S71S7j, i, j = 2, 4, or 6. Two illustrations pertaining to bifurcation of stable steady state S51 to an unstable coexistence steady state are provided in Figure 16. For D = 0.2 day1 and A0 = 1000 mg/L, steady states S4 and S6 are not admissible in single reactor operation (Table 5). Steady state S51 therefore undergoes static bifurcation to steady state S52 at the limit point L (Figures 16a and 16b). For D = 0.3 day1 and A0 = 746 mg/L, steady state S4 is admissible and steady state S6 is not in single reactor operation (Table 5). Steady state S51 undergoes static bifurcation to steady state S54 at the limit point L (Figures 16c and 16d). 4.3. Comparison of Stable Coexistence Steady States. For the two-reactor system, the rate of production of the commensal population, Y, is 0

0

RY ¼ q1e Y1 + q2e Y2

ð22Þ

and the rate of production of biogas, P, is 0

0

RP ¼ ðQ1 V1 + Q2 V2 Þ ¼ εðμ21 Y1 V1 + μ22 Y2 V2 Þ 1535

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

Table 5. Steady State Composition in Single Reactora no. 1 D (day1) 0.29

2 0.3

3 0.3

4 0.2

5 0.29

6 0.25

7 0.3032

A0 (mg/L) 746

746

800

1000

1000

600

775

B (mg/L)

443.77 506.67

C1 (mg/L) 3.062 4.174 C2 (mg/L) 11.207 8.222

506.67 166.53 443.77 278.36 529.66 4.174 8.222

0.891 38.5

3.062 11.21

1.585 21.65

5.121 6.701

A1 (mg/L) 260.79 136.48

239.41 830.31 537.0

302.75 135.89

X1 (mg/L) 1.079

2.68

1.405

0.082

0.5

0.491

2.852

C1 (mg/L) 2.68

6.653

3.488

0.204

1.242

1.222

7.08

Y1 (mg/L) NA

0.059

NA

NA

NA

NA

0.047

Y3 (mg/L) NA

NA

NA

NA

NA

NA

0.009048

A2 (mg/L) 63.713 125.884 76.96 X2 (mg/L) 6.211 2.956 C2 (mg/L) 15.427 7.338

18.613 41.47

5.637 21.09 13.995 52.7

13.41 33.29

38.135 132.733 7.367 2.934 18.337 7.284

Y2 (mg/L) 0.295

0.076

0.235

1.229

0.722

0.399

0.052

Y4 (mg/L) 0.101

NA

0.138

0.337

0.527

NA

0.014

N

17

25

25

25

14

49

25

a

NA = not admissible. N = number of steady states admissible for R < Rc1.

ð23Þ

Y and P, operation at asymmetric steady states S15, S35, S51, and S53 at every R where the asymmetric steady states are admissible. 4.4. Variations in Number of Steady States with R. The maximum number of steady states in a single reactor, seven, is attained in region 7 of the DA0 space (Figure 5), which is rather narrow (Figure 5c) for the kinetic parameters under consideration (Table 6). Operation in region 7 guarantees admissibility of the maximum number of steady states, 49, in a two-reactor system in the immediate neighborhood of R = 0. To examine variation in the number of steady states in the two-reactor system with variation in R, we consider as an illustration operation at D = 0.3032 day1 and A0 = 775 mg/L (Table 7). For single reactor operation, the steady state composition is reported in column 7 of Table 5. On the basis of this composition, one can observe that in the subsystem I space, the following clusters are close to each other: (1) C4 and C5, (2) C7 and C8, (3) C2, C3, C6, and C9. It is evident that Rc1 = R14 = 7.5  103 day1, Rc2 = R15 = 7.757  103 day1, Rc3 = R10 = 5.987  104 day1, Rc4 = R5 = 6.83  105 day1, Rc5 = R12 = 1.072  103 day1, and Rc6 = R1 = 2.373311  106 day1. The locally stable coexistence steady states are (i) S15, S35, S51, S53, and S55 for R < R1, (ii) S35, S51, S53, and S55 for R1 e R < R5, (iii) S35, S51, and S55 for R5 e R < R10, (iv) S51 and S55 for R10 e R < R14, and (v) S55 for R g R14. Clusters C7 and C8 merge at R = R1, clusters C5 and C6 merge at R = R14 and clusters C4 and C9 merge at R = R15. As R is increased, the number of admissible steady states decreases, from 49 for R f 0 to 7 for R > R15, by multiples of two as one crosses a limit point (R = Ri, i = 115, Table 7).

For the profiles of coexistence steady states S35, S51, and S53 in Figures 10, 11, and 16, for every R at which these steady states are admissible, Y01 and Y02 for these are less than Y01 and Y02 for steady state S55, namely, Y2 (Y2 in columns 2, 4, and 5, as appropriate, of Table 5). It is therefore evident from eqs 22 and 23 that, as concerns production of commensal population, Y, and biogas, P, the two-reactor system operation at asymmetric steady states S35, S51, and S53 is inferior to that at symmetric steady state S55. Everywhere in regions 25 and 7 of DA0 space, operation at steady state S55 is superior to, with respect to production of

5. DYNAMICS OF TWO-REACTOR SYSTEM A dynamic operation of the two-reactor system will converge at large to a locally, asymptotically stable steady state Sij belonging to a locally stable cluster and reactor 1 operating at a locally stable steady state Si. Because of one-way interaction between the two reactors, the dynamics of the two-reactor system is decided in order by (i) subsystem I dynamics in reactor 1 (A01, B01, and X01), (ii) subsystem II dynamics in reactor 1 (C01 and Y01), (iii) subsystem I dynamics in reactor 2 (A02, B02, and X02), and (iv) subsystem II

Figure 16. Portraits of (a, c) R and C02 and (b, d) R and Y02 for steady state S51 and the steady state it bifurcates to for D = 0.2 day1 and A0 = 1000 mg/L (parts a and b) and D = 0.3 day1 and A0 = 746 mg/L (parts c and d). Legend ij represents steady state Sij. The steady state composition for single reactor operation is reported in columns 4 and 2 of Table 5. Steady state S51 undergoes static bifurcation at limit points L.

Substitution of the steady-state form of eqs A4 and 5 into the above and further rearrangement leads to 0

0

0

0

RP ¼ ε½k2 ðY1 V1 + Y2 V2 Þ + q1e Y1 + q2e Y2 

1536

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

Table 6. Values of Kinetic Parameters in eqs A1A7 k1

0.004 day1

k2

0.004 day1

K1

160 mg/L

K2

0.82 mg/L

KI R

41.85 mg/L 0.5 L/{mg X 3 day}

β

37.8787 g B/g X

γ

2.45 g C/g X

δ

41.3223 g C/g Y

ε

74.54 {L P}/mg Y

μ10

0.4 day1

μ20

0.4 day1

Table 7. Variation in Number of Steady States with Variation in R for D = 0.3032 day1 and A0 = 775 mg/La Ri  104 i 1

(day1)

clusters

steady state pairs

number of steady states,

merging at Ri

Ni, in (Ri1,Ri)

0.02373311 C7, C8 S12S13, S14S15,

49

S16S17 2 3

0.2419 0.293248

C6 C2

S52S56 S42S46

43 41

4

0.657

C9

S43S47

39

5

0.683

C3

S53S57

37

6

0.914

C6

S72S76

35

7

3.634

C9

S63S67

33

8

5.89

C9

S25S27

31

9

5.98

C6

S34S36

29

5.987 8.93782

C3 C2

S35S37 S62S64

27 25

12 10.72

C3

S73S75

23

13 59.285

C2

S24S26

14 75

C5, C6 S31S32, S51S54,

15 77.574

C4, C9 S21S23, S41S45,

10 11

21 19

S71S74 13

S61S65 a

Number of steady states at R = Ri is (Ni + Ni+1)/2. R0 = 0, N16 = 7.

dynamics in reactor 2 (C02 and Y02) (Figure 17). The dynamics of reactor 1 is discussed in section A.3. 5.1. Dynamics of Reactor 2. 5.1.1. Reactor 1 Operating at Steady State S1. The two-reactor system may eventually operate at one the following steady states: S11, S13 (regions 4, 5, and 7 of DA0 space), S15. Steady states S13 and S15 belong to cluster C8 and hence have the same subsystem I composition, (A02, B02, X02), in reactor 2. The subsystem I space for reactor 2 (A02, B02, and X02) is divided into two regions, the basin of attraction for steady state S11 (cluster C1) and that for steady states in cluster C8 (Table 3 and Figure 17). Convergence to cluster C8 guarantees convergence to coexistence steady state S15 for arbitrary C02 and arbitrary positive Y02 in regions 2 and 3 of DA0 space where steady state S13 is unstable. In regions 4, 5 and 7, where cluster C8 contains two locally stable steady states, S13 and S15, the subsystem II space for reactor 2, (C02,Y02), can be divided into two regions, the basin of attraction for steady state S13 and that for steady state S15, with the separatrix passing through the projection

Figure 17. Flowchart depicting convergence to a locally stable steady state of the two-reactor system. The dynamics of the system is decided in order by (i) subsystem I dynamics in reactor 1, (ii) subsystem II dynamics in reactor 1, (iii) subsystem I dynamics in reactor 2, and (iv) subsystem II dynamics in reactor 2. A0 1(0), B0 1(0), X0 1(0) - initial concentrations of A, B, and X, respectively, in bioreactor 1. C01i, Y01i = concentrations of C and Y, respectively, in bioreactor 1 upon convergence to subsystem I state for group G3. A02i, B02i, X02i = concentrations of A, B, and X, respectively, in bioreactor 2 upon convergence of bioreactor 1 to a stable steady state. C02i, Y02i = concentrations of C and Y, respectively, in bioreactor 2 upon convergence to subsystem I state for cluster C3 or C8 and convergence of bioreactor 1 to a stable steady state.

of steady state S17 in this space. The fate of reactor 2 upon convergence to subsystem I state for cluster C8 is decided by subsystem II state, C02i and Y02i. The two-reactor system will eventually operate at steady state S15 if (C02i,Y02i) lies in the basin of attraction of the steady state (Figure 17). The candidate stable steady states for two-reactor system are S11, S13, and S15, with steady states S13 and S15 being admissible for 0 < R e Rc6. The locally stable steady states then are: (1) 0 < R < Rc6  S11 and S15 (regions 2 and 3), S11, S13, and S15 (regions 4, 5 and 7) and (2) R g Rc6  S11 (globally stable with respect to conditions in reactor 2). 5.1.2. Reactor 1 Operating at Steady State S3. In regions 4, 5, and 7 of DA0 space, the two-reactor system may eventually operate at one of the following steady states: S31, S33, and S35. Steady states S33 and S35 belong to cluster C3 and hence have the same subsystem I composition (A02, B02, X02) in reactor 2. Steady state S31 belongs to cluster C5. The subsystem I space for reactor 2 (A02, B02, and X02) is then divided into two regions, the basin of attraction for steady state S31 (cluster C5) and that for steady states S33 and S35 in cluster C3 (Table 3 and Figure 17). Convergence to cluster C5 guarantees convergence to partial washout steady state S31. In the case of convergence to cluster C3, the subsystem II space for reactor 2, (C02,Y02), can be divided into two regions, the basin of attraction for steady state S33 and that for steady state S35, with the separatrix passing through the projection of steady state S37 in this space. The fate of reactor 2 upon convergence to subsystem I state for cluster C3 is decided by subsystem II state, C02i and Y02i. The two-reactor system will eventually operate at coexistence steady state S35 if (C02i,Y02i) lies in the basin of attraction of the steady state and otherwise lead to washout of Y (Figure 17). The candidate locally stable steady states for two-reactor system are S31, S33 and S35, with steady state S31 being admissible 1537

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research for 0 < R e Rc1 and steady state S35 being admissible for 0 < R e Rc3. The locally stable steady states then are (1) 0 < R < min(Rc1, Rc3)  S31, S33, and S35, (2a) Rc1 e R < Rc3  S33 and S35, (2b) Rc3 e R < Rc1  S31 and S33, and (3) R g max(Rc1,Rc3)  S33 (globally stable with respect to conditions in reactor 2). 5.1.3. Reactor 1 Operating at Steady State S5. The tworeactor system may eventually operate at one of the following coexistence steady states: S51, S53 (regions 4, 5, and 7 of DA0 space), S55. Steady states S53 and S55 belong to cluster C3 and hence have the same subsystem I composition (A02, B02, X02) in reactor 2. The subsystem I space for reactor 2 (A02, B02, and X02) is divided into two regions, the basin of attraction for steady state S51 (cluster C5) and that for steady states S53 and S55 (Table 3 and Figure 17). Convergence to cluster C5 guarantees convergence to steady state S51. In regions 2 and 3 of DA0 space, convergence to cluster C3 guarantees convergence to symmetric steady state S55 as steady state S53 is not admissible. In regions 4, 5, and 7, where there are two locally stable steady states, S53 and S55, in cluster C3, the subsystem II space for reactor 2, (C02,Y02), can be divided into two regions, the basin of attraction for steady state S53 and that for steady state S55, with the separatrix passing through the projection of steady state S57 in this space. The fate of reactor 2 upon convergence to subsystem I state for cluster C3 is decided by subsystem II state, C02i and Y02i. The two-reactor system will eventually operate at steady state S53 if (C02i,Y02i) lies in the basin of attraction of the steady state and otherwise operate eventually at steady state S55 (Figure 17). The candidate locally stable steady states for two-reactor system are S51, S53, and S55, with steady state S51 being admissible for 0 < R e Rc1 and steady state S53 being admissible in regions 4, 5, and 7 of DA0 space for 0 < R e Rc4. In regions 2 and 3, the locally stable steady states are (1) 0 < R < Rc1  S51, S55 and (2) R g Rc1  S55 (globally stable with respect to conditions in reactor 2). In regions 4, 5, and 7, the locally stable steady states are: (1) 0 < R < min(Rc1,Rc4)  S51, S53, and S55, (2a) Rc1 e R < Rc4  S53 and S55, (2b) Rc4 e R < Rc1  S51 and S55, and (3) R g max(Rc1,Rc4)  S55 (globally stable with respect to conditions in reactor 2). Specific numerical results for dynamic behavior of bioreactors in series are not discussed here considering the difficulty associated with graphical presentation of these because dimensions of subsystem I (6), subsystem II (4), and the two-reactor system (10). 5.2. Other Operational Aspects. Although up to 49 steady states are admissible for the two-reactor system, up to nine steady states may be stable, with only up to five steady states permitting stable coexistence. The unstable steady states will not be realized in transient operation of the bioprocess, unless feedback control is used to stabilize one or more of these which permit coexistence. The operation of the two-reactor system is simplified due to one way interaction. Depending on where reactor 1 operates (indicated in parentheses), the steady state operations of interest are (1) S15 (S1), (2) S35 (S3, only regions 4, 5, and 7), (3) S51, S53 (only regions 4, 5, and 7), and S55 (S5). Multiple reactors in series are frequently used for conducting chemical reactions. This configuration provides operational flexibility. Achieving the desired goals, such as reactant conversion and selectivity of the desired product, with a single continuous flow chemical reactor is like putting all eggs in one basket. If the reactor goes out of commission, there can be a substantial loss in productivity unless a standby unit (another reactor of the same size) is available. Use of multiple reactors would provide the same performance as that of a single equivalent reactor when all

ARTICLE

reactors are operational. If some reactors are malfunctioning, these can be put out of commission and bypassed, with the remaining reactors being used until the malfunctioning reactors become operational again. The division of labor provided by multiple reactors makes them attractive over single reactor when operational emergencies arise. For continuous flow biological reactors (continuous cell cultures) with sterile feed, cell washout is an omnipresent steady state. For the bioprocess under consideration here involving host and commensal populations, washout of host and commensal populations (total washout) or washout only of commensal population (partial washout) are omnipresent steady states. Coexistence, implying retention of commensal population, is essential for synthesis of target product P, such as biogas in the present case. In light of this, the series bioreactors offer benefits over and above those offered by series chemical reactors. A single bioreactor of volume (V1 + V2) permits operation at a single coexistence steady state (S5). Besides operation at this steady state (S55), compartmentalization of available volume (V1 + V2), into two interconnected portions of sizes V1 and V2 permits, for certain R, operation at two more coexistence steady states, S15 and S51, in regions 2 and 3 of the DA0 space and at four more coexistence steady states, S15, S35, S51, and S53, in regions 4, 5, and 7. Coexistence is permitted in the entire two-reactor system if it is permitted in reactor 1 and only in reactor 2 otherwise. There is no commensal population and no biogas production in reactor 1 for steady states S15 and S35. The system operation at steady state S15 involves mixed culture in the second reactor, with reactor 1 serving merely as a storage tank. Operation at steady state S35 involves pure culture of host population in the first reactor and mixed culture in the second reactor. For steady state S35, reactor 2 operates on its own as concerns retention of Y and biogas production, unless reactor 1 is rescued via appropriate temporary disturbances in operating conditions, including addition of an inoculum of Y, from steady state S3 to steady state S5. The two-reactor system operation in that case is upgraded from steady state S35 to steady state S55. Additionally, for steady state S15, an inoculum of X must be provided to reactor 1 to upgrade its steady state from S1 to S5 and the two-reactor system operation to steady state S55. Steady states S51 and S53 involve mixed culture in both reactors, with reactor 2 being revived from commensal population washout because of interaction with reactor 1. Any cell species (X or Y here) is protected from washout in a bioreactor whenever it is being continuously introduced into the bioreactor through a feed stream. In the present case, each bioreactor has fresh feed stream (devoid of X and Y) and the second bioreactor an exchange stream from the first bioreactor (Figure 2b). Any such cell species would acquire robustness against washout by virtue of the fact that it can always recover to healthy population levels when conditions are suitable. This is the case for X and Y in reactor 2 when coexistence is maintained in reactor 1. Upgrading of the two-reactor steady state from S53 to S55 can be accomplished via appropriate temporary disturbances in operating conditions for reactor 2, including addition of an inoculum of Y. A switch from steady state S51 to steady state S55 can be accomplished by appropriate temporary disturbances in operating conditions for reactor 2, including addition of inocula of X and Y. No disruption in reactor 1 operation is required in switches from steady states S51 and S53 to steady state S55. Only one coexistence steady state, S5, is locally, asymptotically stable in single bioreactor operation. A significant upset in the 1538

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research operating conditions may cause the reactor to transit to an undesired stable total/partial washout steady state (S1 or S3), forcing termination of bioreactor operation. Upon significant upset in the operating conditions, the two-reactor system may continue to maintain mixed culture owing to transition from one locally stable coexistence steady state to one of the other locally stable coexistence steady states. The risk of bioreactor failure due to permanent loss of Y is therefore significantly reduced. Among the stable coexistence steady states, steady state S55 is the most productive as concerns generation of Y and P. Certain upsets in operating conditions for reactor 2, with the exclusion of the exchange stream between the two reactors, may force the system to transit to a less productive coexistence steady state, S51 or S53. If the upsets are temporary, operation at steady state S55 can be restored by implementing appropriate temporary changes in operating conditions for reactor 2. Certain upsets, at least in the operating conditions for reactor 1, may force the two-reactor system to transit from steady state S55 to a less productive coexistence steady state, S15 or S35. If the upsets are temporary, operation at steady state S55 can be restored by implementing appropriate temporary changes in operating conditions for reactor 1. The two-reactor system is therefore more robust visa-vis single reactor as concerns maintenance of mixed culture.

6. CONCLUSIONS Static and dynamic behavior of a mixed culture in two wellmixed continuous flow bioreactors in series, identical in terms of composition of and dilution rate with respect to fresh feed, has been analyzed in great detail. For the example bioprocess, a single continuous culture may admit up to four coexistence steady states, only one of these being attainable. The series bioreactors may admit up to forty nine steady states, which are comprised of up to forty coexistence steady states, at least at very low interaction rate. The rather large number of steady states arising in the tworeactor system for very low R has been organized into nine clusters, which are divided into symmetric clusters (identical subsystem I composition in two reactors) and asymmetric clusters, and are classified further into three groups. The steady states are divided into symmetric (single reactor) and asymmetric steady states. Every asymmetric steady state (cluster) is not admissible beyond certain R, its limit point, as it will merge with another asymmetric steady state (cluster), with at least one steady state (cluster) in the pair being unstable. The two steady states in each merging pair are of the same type (partial washout or coexistence). The six asymmetric clusters admissible up to a critical R are comprised of three pairs. The two clusters in each pair merge at a limit point. In regions 4, 5, and 7 of the DA0 space, the locally stable cluster C3 is comprised of six asymmetric coexistence steady states, admissible up to a critical R, which are comprised of three pairs. The two steady states in each pair merge at a limit point. Coexistence steady states S15 and S51, which are locally stable and are upgrades of total washout in one bioreactor, undergo static bifurcation to a variety of other unstable coexistence steady states. Only limit point bifurcations are possible for the bioprocess kinetics under consideration. While a single bioreactor can operate at only one coexistence steady state, two bioreactors in series can operate at up to five coexistence steady states over certain ranges of R. The risk of bioreactor failure due to permanent loss of Y is significantly less in a two-reactor system. This system is operationally more flexible and more robust vis-a-vis single reactor

ARTICLE

as concerns maintenance of mixed culture. Emergence of additional steady states demonstrates complexity of the two-reactor system.

’ APPENDIX. BACKGROUND ON SINGLE REACTOR OPERATION For a stand-alone CSTR (Figure 2a), the conservation equations for the species in eq 1 influencing the kinetics are14,34 dA ¼ f1 , f1 ¼ DðA0  AÞ  r1 dt

ðA1Þ

dX ¼ f2 , f2 ¼  DX + ðμ1  k1 ÞX dt

ðA2Þ

dB ¼ f3 , f3 ¼  DB + r1  βμ1 X dt

ðA3Þ

dY ¼ f4 , f4 ¼  DY + ðμ2  k2 ÞY dt

ðA4Þ

dC ¼ f5 , f5 ¼  DC + γμ1 X  δμ2 Y dt

ðA5Þ

with the biogas production rate, Q, being expressed as Q ¼ εμ2 Y

ðA6Þ

The expressions for r1, μ1, and μ2 are34 r1 ¼ RAX, μ1 ¼

μ10 B μ20 C , μ2 ¼ ðA7Þ ðK1 + BÞ ðK2 + CÞð1 + C=KI Þ

The values of kinetic parameters in eqs A1A7 are listed in Table 6.27,34 A.1. Number and Types of Steady States. Three types of steady states are admissible for a single reactor:14 (a) total washout (TW) X = 0, Y = 0; (b) partial washout (PW) X > 0, Y = 0; and (c) coexistence (CO) X > 0, Y > 0. Retention of X at steady state requires that D < (μ10  k1) and A0 g ψðDÞ, ψðDÞ ¼ B +

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 ½βðk1 + DÞ + 4Rβðk1 + DÞB R ðA8Þ

with the expressions for X at steady state (X1 and X2), A, and B, being pffiffiffiffi pffiffiffiffi ϕ Δ ϕ+ Δ , X2 ¼ , X1 ¼ 2Rβðk1 + DÞ 2Rβðk1 + DÞ ϕ ¼ D½RðA0  BÞ  βðk1 + DÞ, Δ ¼ ½ϕ2  4Rβðk1 + DÞD2 B, A¼

DA0 K1 ðk1 + DÞ , B¼ ðD + RXÞ μ10  ðk1 + DÞ

ðA9Þ

It follows from eq A5 that14 (i) C = C [C = γ(k1 + D)X/D] for a partial washout steady state (Y = 0) and (ii) C < C, C = C1, C2, for a coexistence steady state (Y > 0) (Figure 4). When necessary, we will refer to three parameters, Dc (dependent on A0), D*, and A*, 0 which are defined as 

A0 ¼ ψðDc Þ, D ¼ μ2m  k2 , A0 ¼ ψðDÞ, Du ¼ minðDc , DÞ 1539

ðA10Þ

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

We designate steady states for a single reactor using single digit as (notation in parentheses used as an identifier): (S1) TW, (S2) PW1, (S3) PW2, (S4) CO1, (S5) CO2, (S6) CO3, and (S7) CO4 (Table 1, ref 14). On the basis of commonality of subsystem I characteristics, the steady states are divided into three groups: (G1) S1, (G2) S2, S4, S6, and (G3) S3, S5, S7 (Table 1). For the parameters listed in Table 6, coexistence steady states are admissible for 0 < D e D* and A0 g ψ(D) (eq A8). In the DA0 space, the number of physically realizable coexistence steady states will change as one crosses the curves Yk = 0, k = 14.14 Steady states S4 and S6 branch off from steady state S2 and steady states S5 and S7 branch off from steady state S3. The DA0 space is divided into ten regions (Figure 5) depending on which coexistence steady states are admissible (region numbers indicated in parentheses, ref 14): (1) S4, (2) S4, S5, (3) S5, (4) S5, S7, (5) S4, S5, S7, (6) S4, S6, (7) S4, S5, S6, S7, (8) S6, (9) S6, S7, and (10) S7. A.2. Local Stability of Steady States. The local stability of steady states is determined by the eigenvalues of the Jacobian matrix (J) associated with eqs A1A5, J being ( ) ∂fp , p, q ¼ 1  5 , x1 ¼ A, J ¼ jpq , jpq ¼ ∂xq x2 ¼ X, x3 ¼ B, x4 ¼ Y , x5 ¼ C The expression for J is

" J¼

with

2

j11 6 0 A ¼6 4 j31

j12 j22 j32

A D

C B

ðA11Þ

μ1B

ðA12Þ

3 " # 0 j44 j45 7 7 j23 5, B ¼ , C ¼ 0, j54 j55 j33 " # 0 0 0 D¼ ðA13Þ 0 j52 j53

¼ ¼ ¼ ¼ ¼

 ðD + RXÞ, j12 ¼  RA, ðμ1  k1  DÞ, j23 ¼ Xμ1B , j31 ¼ RX, ðRA  βμ1 Þ, j33 ¼  ðD + βXμ1B Þ, ðμ2  k2  DÞ, j45 ¼ Y μ2C , j52 ¼ γμ1 , γXμ1B , j54 ¼  δμ2 , j55 ¼  ðD + δY μ2C Þ, dμ dμ ðA14Þ ¼ 1 , μ2C ¼ 2 dB dC

The characteristic equation for J, |J  λI| = 0 [|M| = det(M) = determinant of matrix M], reduces to jA  λIjjB  λIj ¼ 0

A.3. Dynamics of Reactor 1. The subsystem I space for reactor 1 (A01, B01, and X01) is divided into two regions, one region being the basin of attraction for group 1 (G1) steady state, S1, and the other region being the basin of attraction for group 3 (G3) steady states, S3, S5 and S7 (Table 1, ref 14, Figure 17). The separatrix that separates the two regions passes through the projection of group 2 (G2) steady states, S2, S4, and S6, in subsystem I space. In regions 2 and 3, convergence to subsystem I state for G3 assures convergence to coexistence steady state S5 for arbitrary C02 and arbitrary positive Y02, as the only other admissible steady state in this group, S3, is unstable. In regions 4, 5, and 7, where there are two locally stable steady states, S3 and S5, in group G3, the subsystem II space, (C01,Y01), can be divided into two regions, one region being the basin of attraction for steady state S3 and the other region serving as the basin of attraction for steady state S5. The fate of reactor 1 upon convergence to subsystem I state for group G3 is decided by subsystem II state, C01i and Y01i. Reactor 1 will eventually operate at steady state S5 if (C01i,Y01i) lies in its basin of attraction (Figure 17).

’ AUTHOR INFORMATION

#

The non-zero elements of the Jacobian matrix are j11 j22 j32 j44 j53

and 7. Coexistence steady state S5 is stable in regions 2, 4, 5, and 7 and undergoes Hopf bifurcation at very low A0 and D in region 3 (Figures 5a and 5d). We are not concerned with reactor operation at such low A0 and D and steady state S5 is therefore locally stable in the portion of region 3 of interest.

ðA15Þ

The total washout steady state is locally, asymptotically stable for D > 0. Out of the six steady states involving retention of X (Table 1), only steady states S3 and S5, belonging to group G3, may be stable. The necessary and sufficient condition(s) for local, asymptotic stability, static bifurcation, and Hopf bifurcation of the two steady states are available in ref 14. Partial washout steady state S3 is unstable in regions 2 and 3 and stable in regions 4, 5,

Corresponding Author

*E-mail: [email protected]. Voice: (312) 567-3044. Fax: (312) 567-8874.

’ NOMENCLATURE A0 = concentration of A in the feed to a single reactor and reactor i, i = 1, 2, respectively, mg/L A0*, D*, Dc, Du = defined in eq A10, A0* = ψ(D*), D* = (μ2m  k2), A0 = ψ(Dc), Du = min(Dc,D*) Ai, Bi, Ci, Di = A, B, C, D, respectively, for reactor i, i = 1, 2 AT, BT, CT, DT = defined in eq 19 A, B, C, D = defined in eqs A13 and A14 A, B, C = complex insoluble organics, soluble organics, and volatile acids, respectively C1, C2, Cm = defined in Figure 4, mg/L Cj = steady-state clusters defined in Table 3, j = 19 CO = coexistence steady state X > 0, Y > 0 COj = coexistence steady states with notation in Table 1, j = 14 D = dilution rate for a single reactor, q/V, and reactor i, qi/Vi, day1 Di = dilution rate for reactor i, Di = qi/Vi, i = 1, 2, day1 E = extracellular enzymes synthesized by X fi = defined in eqs A1A5, i = 15 fk2 = defined in eqs 2-6, k = 15 g, x = defined in eq 17 G1, G2, G3 = G1 = TW, G2 and G3 defined in Table 1 J = concentration of specie J, J = A, B, C, X, and Y, in single reactor, mg/L J0 i = concentration of specie J, J = A, B, C, X, and Y, in reactor i, i = 1, 2, mg/L jpq = element of J in row p and column q, p, q = 15, defined in eq A14 when nonzero JT = Jacobian matrix for two coupled reactors defined in eqs 1619 1540

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research J = Jacobian matrix defined in eqs A11A14 k1, k2 = death rate coefficients for acidogens and methanogens, respectively, day1 K1, K2 = saturation constants in growth rate expressions in eq A7 for acidogens and methanogens, respectively, mg/L KI = inhibition coefficient in the growth rate expression for methanogens in eq A7, mg/L P = biogas (metabolites produced by Y) PW = partial washout steady state, X > 0, Y = 0 PWj = partial washout steady states with notation in Table 1, j = 1, 2 q, qi = volumetric feed rate for a single reactor and reactor i, i = 1, 2, respectively, L/day Q, Qi = biogas production rate in a single reactor, defined in eq A6, and reactor i, i = 1, 2, defined in eq 8, L/day q12 = volumetric flow rate from reactor 1 to reactor 2, L/day qie = volumetric effluent rate for reactor i, i = 1, 2, L/day R = interaction rate defined in eq 7, day1 r1 = rate of first step in eq 1, defined in eq A7, mg/{L 3 day} r1i = rate of first step in eq 1, defined in eq A7, in reactor i, i = 1, 2, mg/{L 3 day} Rc1 = R corresponding to limit point for clusters C5 and C6, day1 Rc2 = R corresponding to limit point for clusters C4 and C9, day1 Rc3 = R corresponding to limit point for coexistence steady states S35 and S37, day1 Rc4 = R corresponding to limit point for coexistence steady states S53 and S57, day1 Rc5 = R corresponding to limit point for coexistence steady states S73 and S75, day1 Rc6 = R corresponding to limit point for clusters C7 and C8, day1 Sij = steady state in coupled reactors with reactor 1 operating at steady state Si and reactor 2 operating at steady state Sj when R f 0, i, j = 17 Sj = steady state j in a single reactor, j = 17 t = real time, day TW = total washout steady state, X = Y = 0 V, Vi = culture volume in a single reactor and reactor i, i = 1, 2, respectively, L X1, X2 = defined in eq A9 X, Y = acidogens and methanogens, respectively Yj = Y for COj, j = 14 (Table 1) Greek Symbols

R = kinetic parameter in expression for r1 in eq A7, L/{mg X 3 day} β = yield coefficient for soluble organics utilized by acidogens, g B/{g X} δ = yield coefficient for utilization of volatile acids by methanogens, g C/{g Y} Δ, ϕ = defined in eq A9 ε = yield coefficient for biogas produced by methanogens, {L P}/{mg Y} η = V2/V1, defined in eq 9 γ = yield coefficient for volatile acids produced by acidogens, g C/{g X} λ = eigenvalue of J or JT μ1, μ2 = specific growth rates of acidogens and methanogens, respectively, day1 μ10, μ20 = parameters in eq A7, day1

ARTICLE

μ1B = derivative of μ1 with respect to B μ1i, μ2i = specific growth rates of acidogens and methanogens, respectively, in reactor i, i = 1,2, day1 μ2C = derivative of μ2 with respect to C ψ(D) = defined in eq A8 Subscripts

0 = feed e = reactor effluent i = reactor i, i = 1, 2

’ REFERENCES (1) Birol, I.; Parulekar, S. J.; Teymour, F. Effect of Environmental Partitioning on the Survival and Coexistence of Autocatalytic Replicators. Phys. Rev. E 2002, 66, No. 051916. (2) Chang, S. W.; Baltzis, B. C. Impossibility of Coexistence of Three Pure and Simple Competitors in Configurations of Three Interconnected Chemostats. Biotechnol. Bioeng. 1989, 33, 460–470. (3) Shain, P. Correctly Calculate Conversion in Multichamber Tank Reactors. Chem. Eng. Prog. 1997, 93 (9), 76–80. (4) Katinger, H. W. D.; Scheirer, W. Status and Developments of Animal Cell Technology Using Suspension Culture Techniques. Acta Biotechnol. 1982, 2, 3. (5) Parulekar, S. J.; Hassell, T.; Tripathi, S. C. Recent Developments in Vertebrate Cell Culture Technology. Int. Rev. Cyt. 1992, 142, 145–211. (6) Shuler, M. L.; Kargi, F. Bioprocess Engineering: Basic Concepts, 2nd ed.; Prentice Hall: Upper Saddle River, NJ, 2002. (7) Fredrickson, A. G.; Tsuchiya, H. M. Microbial Kinetics and Dynamics. In Chemical Reactor Theory, A Review; Lapidus, L., Amundson, N. R., Eds.; Prentice Hall: Englewood Cliffs, NJ, 1977; pp 405483. (8) Murray, J. D. Mathematical Biology; Springer-Verlag: Berlin, 1993. (9) Edelstein-Keshet, L. Mathematical Models in Biology; Random House: New York, 1988. (10) Eigen, M.; Schuster, P. The Hypercycle: A Principle of Natural Self-Organization; Springer-Verlag: Berlin, 1979. (11) Parulekar, S. J.; Lim, H. C. Dynamics of Continuous Commensalistic Cultures-I. Multiplicity and Local Stability of Steady States and Bifurcation Analysis. Chem. Eng. Sci. 1986, 41, 2605–2616. (12) Parulekar, S. J.; Lim, H. C. Dynamics of Continuous Commensalistic Cultures-II. Numerical Results for Steady-State Multiplicity Regions and Transient Behavior In-The-Large. Chem. Eng. Sci. 1986, 41, 2617–2628. (13) Zeikus, G., Johnson, E. A., Eds. Mixed Cultures in Biotechnology; McGraw-Hill: New York, 1991. (14) Parulekar, S. J.; Ingle, P. Continuous Commensalistic Cultures with Kinetic Feedback—Static and Dynamic Behavior. AIChE J. 2006, 52, 2949–2963. (15) Sheintuch, M. Dynamics of Commensalistic Systems with SelfAnd Cross-Inhibition. Biotechnol. Bioeng. 1980, 22 (12), 2557–2577. (16) Stephanopoulos, G. The Dynamics of Commensalism. Biotechnol. Bioeng. 1981, 23, 2243–2255. (17) Lee, J.; Parulekar, S. J. Periodic Operation of Continuous Recombinant Cultures Improves Antibiotic Selection. Chem. Eng. Sci. 1996, 51, 217–231. (18) Abashar, M. E. E.; Judd, M. R. Synchronization of Chaotic Nonlinear Oscillators: Study of Two Coupled CSTRs. Chem. Eng. Sci. 1998, 53 (21), 3741–3750. (19) Chen, C.-C.; Fu, C.-C.; Tsai, C.-H. Stabilized Chaotic Dynamics of Coupled Nonisothermal CSTRs. Chem. Eng. Sci. 1996, 51 (23), 5159–5169. (20) Taylor, M. A.; Kevrekidis, I. G. Couple, Double, Toil and Trouble: A Computer Assisted Study of Two Coupled CSTRs. Chem. Eng. Sci. 1993, 48, 2129–2149. (21) Baltzis, B. C.; Fredrickson, A. G. Competition of Two Microbial Populations for a Single Resource in a Chemostat When One of Them Exhibits Wall Attachment. Biotechnol. Bioeng. 1983, 25, 2419–2439. 1541

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542

Industrial & Engineering Chemistry Research

ARTICLE

(22) Kung, C.; Baltzis, B. C. Operating Parameters’ Effects on the Outcome of Pure and Simple Competition Between Two Populations in Configurations of Two Interconnected Chemostats. Biotechnol. Bioeng. 1987, 30 (9), 1006–1018. (23) Pavlou, S.; Kevrekidis, I. G.; Lyberatos, G. On the Coexistence of Competing Microbial Species in a Chemostat Under Cycling. Biotechnol. Bioeng. 1990, 35, 224–232. (24) Stephanopoulos, G.; Fredrickson, A. G. Effect of Spatial Inhomogeneities on the Coexistence of Competing Microbial Populations. Biotechnol. Bioeng. 1979, 21, 1491–1498. (25) Reith, J. H.; Wijffels, R. H.; Barten, H., Eds. Bio-Methane & BioHydrogen: Status and Perspectives of Biological Methane and Hydrogen Production; Dutch Biological Hydrogen Foundation, 2003. (26) Simeonov, I.; Momchev, V.; Grancharov, D. Dynamic Modeling of Mesophilic Anaerobic Digestion of Animal Waste. Water Res. 1996, 30, 1087–1094. (27) Simeonov, I.; Stoyanov, S. Modeling and Dynamic Compensator Control of the Anaerobic Digestion of Organic Wastes. Chem. Biochem. Eng. Quart. 2003, 17 (4), 285–292. (28) Hobson, P. N.; Wheatley, A. D. Anaerobic Digestion: Modern Theory and Practice. Elsevier Applied Science: London, 1992. (29) Wheatley, A. D., Ed. Anaerobic Digestion: A Waste Treatment Technology.: Elsevier Applied Science: London, 1991. (30) Bouskova, A.; Dohanyos, M.; Schmidt, J. E.; Angelidaki, I. Strategies for Changing Temperature from Mesophilic to Thermophilic Conditions in Anaerobic CSTR Reactors Treating Sewage Sludge. Water Res. 2005, 39, 1481–1488. (31) Dinsdale, R. M.; Premier, G. C.; Hawkes, F. R.; Hawkes, D. L. Two-Stage Anaerobic Co-Digestion of Waste Activated Sludge and Fruit/Vegetable Waste Using Inclined Tubular Digesters. Bioresour. Technol. 2000, 72, 159–168. (32) Hill, D. T.; Barth, C. L. A Dynamic Model for Simulation of Animal Waste Digestion. J. Wat. Poll. Control Fed. 1977, 10, 2129–2143. (33) Keshtkar, A.; Meyssami, B.; Abolhamd, G.; Ghaforian, H.; Khalagi-Asadi, M. Mathematical Modeling of Non-Ideal Mixing Continuous Flow Reactors for Anaerobic Digestion of Cattle Manure. Bioresour. Technol. 2003, 87, 113–124. (34) Noykova, N.; Muller, T. G.; Gyllenberg, M.; Timmer, J. Quantitative Analyses of Anaerobic Wastewater Treatment Processes: Identifiability and Parameter Estimation. Biotechnol. Bioeng. 2002, 78 (1), 89–103. (35) Pahl, M.; Lunze, J. Dynamic Modelling of a Biogas Tower Reactor. Chem. Eng. Sci. 1997, 53, 995–1007. (36) Tartakovsky, B.; Guiot, S. R.; Sheintuch, M. Modeling and Analysis of Co-Immobilized Aerobic/Anaerobic Mixed Cultures. Biotechnol. Prog. 1998, 14, 672–679. (37) Jeyaseelan, S. A Simple Mathematical Model for Anaerobic Digestion Process. Water Sci. Technol. 1997, 35, 185–191. (38) Angelidaki, I.; Ellagard, L.; Ahring, B. K. A Comprehensive Model of Anaerobic Bioconversion of Complex Substrates to Biogas. Biotechnol. Bioeng. 1999, 63, 363–372. (39) Husain, A. Mathematical Models of the Kinetics of Anaerobic Digestion—A Selected Review. Biomass Bioenergy 1998, 14, 561–571. (40) v. Munch, E.; Keller, J.; Lant, P.; Newell, R. Mathematical Modelling of Prefermenters-I. Model Development and Verification. Water Res. 1999, 33, 2757–2768. (41) Dhooge, A.; Govaerts, W.; Kuznetsov, Y. A.; Mestrom, W.; Riet, A. M. CL_MATCONT: a Continuation Toolbox in Matlab. In SAC ’03: Proceedings of the 2003 ACM Symposium on Applied Computing; ACM Press: Melbourne, FL, 2003; pp 161166.

1542

dx.doi.org/10.1021/ie200103g |Ind. Eng. Chem. Res. 2012, 51, 1525–1542