Competition between Robust and Nonrobust Autocatalytic Replicators

3630

Ind. Eng. Chem. Res. 2002, 41, 3630-3641

Competition between Robust and Nonrobust Autocatalytic Replicators Woraritt Chaivorapoj, I4 nanc¸ Birol, and Fouad Teymour* Department of Chemical and Environmental Engineering, Illinois Institute of Technology, Room 127, 10 W. 33rd Street, Chicago, Illinois 60616

The competition of two species for a common resource is illustrated using the paradigm of autocatalytic replicators inhabiting a continuous stirred tank reactor (CSTR) environment that is continuously fed with the resource. In most cases presented, one species is robust (appears in reactor feed) while the other is not. The introduction of the second (invading) species into the CSTR via an unsustained disturbance has a strong effect on the steady-state and dynamic behavior of the first (host) species. New steady states are added to the bifurcation diagram that show that the invading species can coexist in the system with the host species when its growth and death characteristics are similar to those of the latter. The population levels of the host species are greatly reduced in these cases as a result of the considerable decrease in resource concentration at steady state. Open-loop strategies for the elimination of the invading species are developed and discussed. These strategies involve the manipulation of the reactor residence time to destabilize the states of coexistence. 1. Introduction The isothermal continuous stirred tank reactor (CSTR) with an autocatalytic reaction has been the subject of interest for over two decades,1-14 mostly because it combines mathematical simplicity with a richness of steady-state and dynamic behavior. The complexity of the system is the result of the autocatalytic feedback, where the concentration of the product of the reaction affects its production rate. The earliest investigations of general autocatalytic reactions were performed by Lin.1,2 In a theoretical study,1 Lin analyzed autocatalytic reactions of type A + B f 2B in a CSTR for steady-state multiplicity and stability, showing that the parameters of importance to multiplicity criteria are the reaction orders, the dimensionless residence time, and the ratio of feed concentration of B to that of A. He showed that multiplicity is guaranteed by limiting the dimensionless residence time in a proper range, in addition to ensuring that necessary conditions for the reaction orders and the ratio of the feed concentrations are fulfilled. However, the derived multiplicity criterion can be applied only to the case when there is some amount of autocatalytic species present in the feed. A second article2 completed the framework by establishing ranges of multiplicity, stability, and dynamics for a generic autocatalytic reaction mA + nB f (n + 1)B in an isothermal CSTR where there is no autocatalyst in the feed. It is well-known that the initial conditions will influence the yield or the selectivity of the desired product in the presence of a competing reaction in the CSTR.3,4 Heinrichs and Schneider4 adapted Schlo¨gl’s analysis of the CSTR and demonstrated the phenomenon of “slowing down” where the relaxation times for the decay of small perturbations could become very long compared to the residence time of the reactor. Because * Corresponding author. Tel.: +1-312-567-8947. Fax: +1312-567-8874. E-mail: [email protected].

neither of these works includes the decay reaction (B f inert) in the system, the concentrations of A and B are not independent, but rather are uniquely related to the inlet composition by the net reaction stoichiometry (A f B). Consequently, this one-variable system could not display any complex behavior such as oscillatory solutions. The works of Gray and Scott on quadratic autocatalysis (A + B f 2B)7 and cubic autocatalysis (A + 2B f 3B)8 with decay (B f inert) in an isothermal CSTR, on the other hand, show exotic behavior, such as isolas and mushrooms in plots of conversion (γ) versus residence time (τ). Furthermore, in these systems, the analytical analysis of stability is still tractable. For cubic autocatalysis with no autocatalyst B in the inflow, a maximum of three steady states is possible. On the (γ, τ) plane, the range of multiple solutions is found on a closed curve (or isola) that lies above the trivial solution that exists for all residence times and corresponds to zero conversion of A to B. The trivial solution is globally stable, and small perturbations around it decay monotonically. As some autocatalyst is added to the feed in increasing concentration, the trivial solution is distorted upward until it eventually coalesces with the bottom of the isola, yielding a mushroomshaped curve. The middle steady-state solution on both sides of this mushroom is always unstable and is a saddle point. The high-conversion steady-state solution has ranges of stability, where perturbations decay via damped oscillations, but it also becomes unstable after a Hopf bifurcation point is crossed. On many occasions, one can draw analogies between reaction kinetics and population dynamics.15 For example, quadratic autocatalysis is analogous to the reproduction of single cells by division,16-20 and cubic autocatalysis is similar to sexual reproduction in higher organisms.19,20 It should be noted that the population dynamics analogy of autocatalytic reactions could also be carried over to DNA replication, RNA polymeriza-

10.1021/ie010302p CCC: $22.00 © 2002 American Chemical Society Published on Web 06/28/2002

Ind. Eng. Chem. Res., Vol. 41, No. 15, 2002 3631

tion,21-23 and other enzymatic reactions in biosystems14 or to stock market dynamics.24-26 However, to the best of our knowledge, no extensive study has been performed on multiple autocatalators competing for the same resource, i.e., in pure and simple competition as defined by Fredrickson and Stephanopoulos28 for mixed microbial cultures. It is shown both analytically using Monod kinetics and experimentally that, in a homogeneous environment with time-invariant input, two microbial populations competing for a single resource cannot coexist in a steady state27,28 if neither is fed into the system. However, coexistence of the competing microbials is possible if an inhomogeneity is introduced into the environment (e.g., by wall attachment29). In a study of multiple cubic autocatalysis, in agreement with the Monod kinetics results, Birol and Teymour30 found that, in the competition of multiple autocatalysts for a shared resource in an isothermal CSTR, stable steadystate coexistence is not possible with time-invariant input and no autocatalysts in the feed. It is evident then that the long-term survival of any species in this type of CSTR environment is guaranteed if this specific species is present in the reactor feed. Even small concentrations of the autocatalytic species in the feed will ensure that, whenever the conditions in the environment become favorable, the species will thrive in the system. For the purpose of this study we will associate this circumstance with the trait of robustness. The species thus is considered “robust” in as much as it can never be “totally” removed from the system. In this work, we study a homogeneous isothermal CSTR with two cubic autocatalytic reactions, where both autocatalysts exhibit decay and where at least one of the autocatalysts is robust. When only one robust species is considered, it will be termed the “host species”, and the other autocatalyst will be considered an “invading species”. This second species will only be introduced to the reactor by means of a disturbance in the feed. Unlike nonautocatalytic reactions, the introduction of a second species leads to permanent changes in the static structure of the system, in addition to altering the dynamic response. The purpose of this work is to gain an understanding of the statics and dynamics of competition of two cubic autocatalytic reactions in an isothermal CSTR with one robust species in the environment. In the next section, we present the reaction mechanism for the case of two competing cubic autocatalysts in an isothermal CSTR and introduce the mathematical model equations. In the following section, we first summarize the work of Gray and Scott8 for the host species and then investigate the altered bifurcation diagram of the system after the introduction of the invading species in the absence and presence of the latter in the feed. In section 4, we present case studies for the dynamic response of the system and propose open-loop control strategies to return the system to a “healthy” operating point, i.e., to remove the invading species from the reactor.

krp

R + 2P 98 3P kdp

P 98 Dp krq

R + 2Q 98 3Q kdq

Q 98 Dq

(1) (2) (3) (4)

where Dp and Dq are the dead forms of species P and Q, respectively, and krp (krq) and kdp (kdq) are the kinetic rate constants for the reproduction and death reactions of species P (Q), respectively. Because the second species Q utilizes the same resource R as the first species P, the interaction between the two is considered pure and simple competition. In general, the resource R (at a concentration Ro) and a trace amount of species P, with concentration Po, are fed continuously to the reactor. As Q is introduced into the reactor as part of a transient pulse disturbance in the feed, we will consider the concentration of the second species in the feed, Qo, to be equal to zero, unless stated otherwise. Note that the static structure of the system is altered by this transient disturbance, as we will demonstrate later. The model equations are written in terms of the dimensionless concentrations r ) R/Ro, p ) P/Ro, po ) Po/Ro, and q ) Q/Ro; the dilated time t ) krpRo2t′; the scaled parameters dp ) kdp/(krpRo2), dq ) kdq/(krpRo2), and kq ) krq/krp; and the dilated residence time T ) krpθRo2. These equations are given by

1 dr ) -rp2 - kqrq2 + (1 - r) dt T

(5)

1 dp ) rp2 - dpp + (po - p) dt T

(6)

dq 1 ) kqrq2 - dqq - q dt T

(7)

At steady state, where dr/dt ) dp/dt ) dq/dt ) 0, there will be zero and nonzero steady-state solutions of q in eq 7. The system will represent single-species autocatalysis for q ) 0 and two-species autocatalysis for q * 0. 3. Steady-State Analysis 3.1. Single-Species Autocatalysis. We first consider the system with no invading species Q (q ) 0) in the reactor. Thus, the dimensionless mass balance equations reduce to the system investigated by Gray and Scott,8 namely

1 dr ) -rp2 + (1 - r) dt T

(8)

dp 1 ) rp2 - dpp + (po - p) dt T

(9)

2. Mathematical Model The system considered is an isothermal CSTR in which a resource, R, is utilized by two distinct cubic autocatalytic species P and Q for their reproduction. Thus, the reaction scheme in the CSTR is given by

Their analysis proposed a two-parameter space diagram in which po and dp were used to map different types of bifurcation diagrams to different regions in the parameter space, with T being used as the primary bifurcation parameter. Figure 1 shows the partitioning of this

3632

Ind. Eng. Chem. Res., Vol. 41, No. 15, 2002

T2kq2r4 - 2T2kq2(po + 1)r3 + Tkq[kq(1 + dpT)2 + 2(1 + dqT)2 + Tkq(po + 1)2]r2 - Tkq[kq(1 + dpT)2 + 2(po + 1)(1 + dqT)2]r + (1 + dqT)2[kq(1 + dpT)2 + (1 + dqT)2] ) 0

Figure 1. (dp, po) parameter space is partitioned into three regions by curves 1, 2, and 3. Regions M, I, and S define a mushroom, an isola, and a single solution in the bifurcation diagram, respectively. Curves 1, 2, and 3 meet at a cusp point with coordinates po ) 1/8, dp ) 33/44. Points a, b, and c are representative of the three types of bifurcation diagrams as illustrated later.

(po, dp) parameter space, presenting three regions M, I, and S in which different types of bifurcation diagrams arise. Regions M, I, and S define a mushroom, an isola, and a unique solution in the bifurcation diagram, respectively. For high concentrations of the species in the feed (po > 1/8), there will be only one real steadystate solution for eqs 8 and 9 for all dimensionless residence times, regardless of the dimensionless death rate of this species (dp). For lower po, either one or three real solutions are possible, depending on the value of dp. Curves 1 and 3 in Figure 1 separate the region of single solutions (S) from the regions of multiple solutions (M and I), whereas curve 2 indicates the transformation from an isola (region I) to a mushroom (region M). All of these curves meet at a cusp point at po ) 1/8, dp ) 33/44, which represents a higher-order singularity. 3.2. Multiple-Species Autocatalysis. Now, consider the two-species autocatalysis (q * 0), where the second species is introduced into the single-species system as a pulse disturbance in the feed. From eq 7, the nonzero steady-state dimensionless concentration q can be found as

1 + dq T q) kq r

Note that the coefficients of the polynomial have alternating signs, which is a necessary condition for the steady-state R concentrations to be positive. Now, for any given set of system parameters kq, T, dp, dq, and po, there will be either zero, two, or four real roots for eq 12. Thus, the system can have a maximum of four steady states for the two-species case (q * 0), in addition to the maximum of three steady states resulting from the single-species case (q ) 0). We can further generalize this analysis to a system with a robust species, p, and n invading species qi, i ) 1, ..., n, of corresponding concentration qi. The nontrivial steady-state concentration of each invading species is given by

1 + dqi T qi ) kqir

(13)

and the corresponding robust-species concentration can be found as

p)

x

-

∑i

(

1

T

)

+ dqi

2

+

rT

r4 - 2(1 + po)r3

[

(10)

x

kqr2

12 2 r T

)] 1 - [S(1 + p ) + T(d + ) ]r T 1 + S[S + T(d + ) ] ) 0 T (

+ S + (1 + po)2 + T dp +

2

p

2

p

p)

(14)

Then, the steady-state resource concentrations still come from the roots of a fourth-order polynomial

Substituting eq 10 into eq 5 yields the dimensionless concentration p as

q

1-r

kqir2

o

1 +d ) ( T -

(12)

(15)

2

+

(1 - r) rT

where

(11)

Note that the r terms in the denominators of eqs 10 and 11 will not impose any singularity, as the steady-state value of r cannot be 0, because a constant supply of r is fed to the system. By substituting eq 11 into eq 6, the steady-state dimensionless concentrations r are given by the roots of the fourth-order polynomial

∑i

S)T

(

1

T

)

+ dqi kq i

2

(16)

Note that, when n ) 0, the polynomial order reduces to three, because r ) 0 is not a solution. Thus, a system with n invading species has 2n possible subsystems with


Figure 2. Bifurcation diagrams in the (T, r), (T, p), and (T, q) planes of the two-species system with constant po ) 0.04, dp ) 0.03, and kq/dq ) 64 show the effect of changing dq on the location of the isola in the (T, r) and (T, q) planes. dq is (a) 6.1440, (b) 0.7680, (c) 0.3840, (d) 0.1920, (e) 0.0480, (f) 0.0240, (g) 0.0030, (h) 0.0015, and (i) 1.875 × 10-4. Note that the lines for zero solution of Q are not shown in the (T, q) plane.

zero, one, two, ..., and n nonzero invading-species concentration values. Each of these subsystems, except the first, has up to four distinct solutions, and the system with no invading species has up to three distinct steady states. Therefore, a system with one robust species and n invading species has up to 2n+2 - 1 steady states. In this study, as in the work of Gray and Scott,8 the dimensionless residence time (T) is selected as the main bifurcation parameter. Because of the difficulties of solving the 4th-order parametric polynomial of eq 12 analytically, the CONTENT software31 is used to compute the bifurcation diagrams of the resource R, species P, and species Q concentrations versus T. This software makes use of bifurcation methods along with numerical continuation techniques in tracing branches of steady state, as well as periodic solutions of sets of ordinary differential equations. For example, in Figure 2a, the bifurcation diagram of the two-species autocatalysis is shown for the parameter set po ) 0.04, dp ) 0.03, dq ) 6.1440, and kq ) 393.216 (64dq), where we see the

emergence of new steady states, characterized by an isola in the (T, r) plane. The solid curves in Figure 2a show the stable steady states, and the dashed curves the unstable ones. As usually encountered in multiplicity regions, branches of stable nodes are separated by a branch of unstable saddles. Additionally, a dynamic exchange of stability occurs at Hopf bifurcation points (seen at the transition from stability to instability on the nonsaddle branches) at which a stable focus gives way to an unstable focus leading to a (stable or unstable) limit cycle. The transition from a stable node to a stable focus occurs, as expected, in the near vicinity of the Hopf points before these points are crossed. It can also be seen that the interaction isola (Figure 2c and g) is fully unstable and consists of a saddle branch and an unstable nodal branch. Note that the mushroom curves in the (T, r) and (T, p) planes represent the single-host-species case (i.e., q ) 0). Because of the competition for the single resource R, with species Q, the reproduction performance of species P in the two-species case is decreased consider-

3634


ably, so that its isola in the (T, p) plane is flattened near low concentration levels and is hardly seen in Figure 2a, because it is trapped beneath the lower residence time branch of the mushroom. This is because the reproduction performance of species P in the two-species case cannot be better than that in the single-species case. The isolas in the (T, r) and (T, q) planes lie in the same range of residence times on the bifurcation diagrams as they are different projections of the same isola. Furthermore, species Q will not be able to exist in the system if the residence time of the system is not in the range for which its isola exists. The effects of changing the system parameters dq and po are investigated next. Effect of dq. The effect of decreasing dq while keeping the ratio of kq/dq and other system parameters constant is illustrated in Figure 2. Mushroom curves of the single-species steady state in the (T, r) and (T, p) planes in Figure 2a-i correspond to point a in Figure 1, where po ) 0.04 and dp ) 0.03. It should be noted that the size of the isola in Figure 2a is primarily determined by the ratio of kq to dq,30 which is fixed at 64 in this study. Keeping the ratio of kq to dq constant, the isola of Figure 2a is translated to a higher residence time in the (T, r) plane by decreasing the value of dq from 6.1440 to 0.7680 (Figure 2b). As dq is decreased further to 0.3840, the isola moves further inside the mushroom of the host species, and we observe the birth of a new isola in the bifurcation diagram, because of the interaction of two species. Note that there are no stable steady states on this interaction isola, just as in the variant of the problem with po ) 0.30 For dq between 0.1920 and 0.0240 (Figure 2d and f), the interaction isola disappears from the bifurcation diagram, as the two isolas merge onto a single distorted isola. Then, around dq ) 0.0030, the interaction isola separates again (Figure 2g), and around dq ) 0.0015, it vanishes (Figure 2h). After passing through the mushroom, the complete isola of the two-species steady states is seen again at the other side of the mushroom (Figure 2i). Notice that, much like the single-species steady states, the stable region for two-species steady states is always on the lower part of the isola on the (T, r) plane. From Figure 2a-i, it is obviously seen that the major effect of changing dq while keeping kq/dq constant is on the location of the isola and that the observed deformations are due to the interaction between the two species. From the host species’ point of view, it is desirable to operate the system with a residence time in the range where the concentration p is ignited on the higher part of the mushroom. For the mushroom of the singlespecies curve in Figure 2a, this range is approximately 5-120. We will call this the “range of survival”. The residence time ranges of the isolas of the two-species steady states in Figure 2a (dq ) 6.1440) or 2i (dq ) 0.0001875) are not in this range because the death rate (dq) is much higher or the reproduction rate (kq) is much lower than those of species P, respectively. This implies that the flow rate in the range of survival is too low for high-dq species and too high for low-dq species to survive. These kinds of invading species will eventually be starved to death for high dq and washed out for low dq. However, they can alter the steady state of the system long after their disappearance whenever multiple steady states exist in the parameter space of the host species. An example for this case is presented in section 4. On the other hand, the invading species Q

Figure 3. Two-parameter continuation diagram of limit points for the two-species system in the (dq, T) parameter space. Other parameters are set at the same values as were used in Figure 2.

will survive in the range of operating conditions desirable for the host species when its death rate is comparable to that of P, as in the cases of parts b-h of Figure 2. Coexistence of the two species results in the intersection of the residence time ranges of the stable region of the isola and the host-species range of survival. The population level of the host species P is therefore decreased dramatically by the presence of these kinds of invading species Q. Parameter Space Classification. Figure 2 illustrates the transition between different bifurcation diagram (BD) varieties as the dq parameter is varied from high to low values. The information provided by this figure about the critical transition values could be portrayed in a more compact manner on a twoparameter continuation diagram, as in Figure 3. This figure depicts the variation of the loci of all limit points (LPs) in the primary bifurcation parameter T as a result of the variation of the secondary bifurcation parameter dq. In Figure 3, lines 1 and 2 correspond to the outermost limit points on the large two-species isola, whereas closed curves 3 and 4 represent the appearance/ disappearance of the smaller interaction isola and the joining/disjoining of the two isolas into a “double-limit” isola, respectively. Vertical lines 5-8 correspond to the four limit points on the host-species mushroom. As expected, these points are not affected by any of the properties of the invading species. When scanning the diagram from top to bottom, one would encounter regions where the bifurcation diagram displays 6, 8, 10, 8, and 6 limit points as dq is lowered. Each of these transitions marks a qualitative change in the shape of the bifurcation diagram and thus provides means for classifying variety. Further transitions can be seen when the LP loci intersect the vertical loci 5-8 corresponding to the mushroom curve. Additional classification could also be obtained by the two-parameter continuation of Hopf bifurcation points, but this approach was not attempted here. As the values of parameters other than the two bifurcation parameters change, not only would the transitions occur at different values of dq, but also other types of behavior could be encountered. To investigate this, we constructed parameter space classification diagrams as in Figures 4 and 5. There, a third parameter is allowed to vary to locate various regions of bifurcation variety in a two-dimensional parameter space. The parameter space selected consists of the two


Figure 4. Parameter space classification diagram for the twospecies system in (dq, kq/dq) space. Each region reflects a different variety of bifurcation in the main bifurcation parameter T. Figure 6. Typical bifurcation diagram in region IV-f of the parameter space classification (Figure 5b). kq/dq ) 33, dq ) dp ) 0.03 and other parameters are as in Figure 2. Steady-state stability is not indicated on this figure.

Figure 5. (a) Further classification of parameter space according to the location of isola LPs relative to the LPs of the mushroom of the host species. (b) Enlarged view at lower values of kq/dq shows additional regions that were not observed under the conditions of Figure 2.

bifurcation parameters kq/dq and dq, and the bifurcation variety is described with respect to the main free parameter T. Note that we have elected to vary the ratio kq/dq rather than the individual value of kq. This choice was guided by the findings of Birol and Teymour,30 who concluded (for a system of nonrobust species) that changing dq while maintaining kq/dq constant preserves the size of the isola on the logarithmic residence time scale and merely translates it to higher values as dq is lowered. This conclusion obviously carries over for the

system under investigation here as illustrated by Figures 2 and 3. If, however, the value of po is increased significantly, one would expect substantial deviations from this behavior. The parameter space classification in Figure 4 results when limit points on isolas are the only ones traced. Four distinct regions are encountered, which either have no isolas, a single isola, two distinct isolas, or a single double-limit isola. It is interesting to note that the region of highest activity of interaction between the two species is centered around the point at which their death rates are equal (dq ) dp ) 0.03) and that the diagram is symmetric about this value (on a logarithmic dq scale). It is also interesting to note that the condition for steady-state survival of the invading species in the system mimics closely that derived by Birol and Teymour,30 which requires that kq/dq be larger than 16 for any nonrobust species Q to survive in the system. A close examination of Figure 3 reveals slight deviations from this condition around the plane of symmetry, where a maximum on the “almost” horizontal line separating regions I and II is observed. This is caused by the robustness of the competing species, which makes it tougher for the invader to coexist (i.e., requires a slightly higher critical limit for kq/dq). A further refinement of the parameter space is presented in Figure 5. This was achieved by tracking the loci of the intersections between lines 1 and 2 of Figure 3 and the LPs of the mushroom curve. The bifurcation diagram in each subregion shown in Figure 5 essentially differs from the others in terms of the relative location of its isolas in relation to the mushroom. For example, the bifurcation diagrams of Figure 2a-i fall in regions II-a, II-c, III-a, IV-a, IV-c, IV-c, IIIe, II-d, and II-f, respectively. In all of these cases, the outer isola size is constant (because kq/dq is constant), and it is larger than the width of the mushroom curve. As kq/dq is decreased further, the outer isola shrinks to sizes that would completely fit inside the mushroom or even fully under one of its eaves. This gives rise to additional subregions that cannot be observed under the conditions of Figure 2. These regions are shown in the enlarged parameter space section of Figure 5b. A typical bifurcation diagram in one of these new regions is shown

3636


Figure 7. Bifurcation diagrams in the (T, r), (T, p), and (T, q) planes of the two-species system with constant dp ) 0.08, dq ) 0.0240, and kq/dq ) 64 show the effect of changing po on the structure of the steady state of the system. po is (a) 0.04, (b) 0.08, (c) 0.12, and (d) 0.16.

in Figure 6. In this diagram, full symmetry (in log T) is obtained because dq and dp are equal. It is interesting to observe that, in this figure (and any others discussed previously), the inner limit points on the double-limit isola are always confined inside the mushroom gap. This behavior is confirmed by line 4 of Figure 3 and is seen generically throughout. Effect of po. To observe the effect of the amount of the host species present in the feed, po, we allow its concentration to vary from 0.04 to 0.16 while keeping the other system parameters constant at dp ) 0.08, dq ) 0.0240, and kq ) 1.536 (or 64dq). These results are illustrated in Figure 7 and show that the interaction between the two species develops as po is increased. When present at relatively low levels (po ) 0.04, Figure 7a) in the feed, the host species has a unique steady state, reflecting a low concentration in the system. Naturally, whenever the invading species Q is introduced, it dominates over the range of residence times defined by the stable portion of its isola. Notice the thin isola tucked underneath the single steady-state curve on the (T, p) bifurcation diagram of Figure 7a, reflecting a decreased concentration for p and corresponding to a high concentration of q as in the (T, q) diagram. As po is increased slightly (Figure 7b), the host species can better survive in the system, as reflected by the appear-

ance of its own isola with higher concentration levels. The isola of the invading species is almost unaffected at this stage. Yet, as the feed concentration of p is further increased, as in Figure 7c,d, the invading-species isola is severely distorted. As a mushroom-shaped curve develops for p (Figure 7c), the underside of the q isola seems to echo the curvature of the mushroom. Even more distortion is seen after the limit points of the mushroom curve disappear, and the host-species space reverts to a single steady state (Figure 7d). It should be noted, however, that all of these changes have a minimal direct impact on the survival fate of the invading species, as the range and location of its stable branch on the isola are almost unaffected. On the other hand, the fate of species Q is indirectly impacted by the increased survival ability of the host species, as is evident from the (T, p) bifurcation diagrams. It becomes plausible then that, in a system infected by Q, the manipulation of po can lead to a reversal of the dominance of the species competing for survival in the environment, thus enabling the removal of the invading species. It should be noted that, in the case of Figure 7, the reproduction performance of the invading species Q was selected to be superior to that of the host species P, as kq/dq > kp/dp. Thus, even if the invading species Q is


Figure 8. Bifurcation diagrams of the two-robust-species system in the (T, r), (T, p), and (T, q) planes for the system parameters in Figure 2c are shown. The new system parameter, qo, is (a) 0.01, (b) 0.04, and (c) 0.1. The solid lines represent the stable regions, and the dashed lines the unstable ones.

allowed to dominate whenever it is introduced into the system, an increase in the feed concentration of species P (po) can lead the host species P to dominate the system instead. A consideration of the trends in Figure 7a-d is essential and instrumental whenever an attempt is made to devise a strategy for elimination of the invasion. 3.3. Multiple-Robust-Species Autocatalysis. In all of the cases considered so far, robustness was bestowed only on the host species, ensuring that the elimination of the invading species is possible. When two competing species are both robust, permanent coexistence in the system prevails over a range of conditions, and interesting interactions result. To analyze these phenomena, we briefly investigate the system in which species Q is also fed continuously to the reactor. To include this effect in the model, we have to modify eq 7 to read

1 dq ) kqrq2 - dqq - (qo - q) dt T

(17)

where qo ) Qo/Ro is the dimensionless feed concentration of species Q. Now, total washout (q ) 0) is no longer one of the steady-state solutions of Q, but a maximum of seven steady states is still possible for this two-robustspecies system, as shown in Figure 8. To illustrate the effect of changing qo on the steady states of the system, the same set of system parameters as in Figure 2c [po ) 0.04, dp ) 0.03, dq ) 0.3840, and kq ) 24.576 (64dq)]

was used. Three values of qo are used (0.01, 0.04, and 0.1, corresponding to values less than, equal to, and greater than po, respectively) to obtain the bifurcation diagrams in Figure 8. At qo ) 0.01 (Figure 8a), the bifurcation diagrams in the (T, r), (T, p), and (T, q) planes are similar to those in Figure 2c, indicating that the level of qo is too low to alter the steady-state structure of the system significantly. For qo ) 0.04 (Figure 8b), the level of qo is high enough to open the isola of species Q to start a mushroom at lower T, but its higher-T branch merges with the lower-T branch of the mushroom of species P, resulting in an interaction double mushroom in the (T, r) plane. In the (T, p) and (T, q) planes, the isola in Figure 8a opens to distort the original mushroom. As qo is increased to 0.1, the distorted mushroom in the (T, r) plane widens (Figure 8c), and the residence time range seems to be divided into three distinct regions where either q dominates, the two species coexist at comparable levels, or p dominates as T is increased. It is obvious from Figure 8a-c that the effect of increasing qo is similar to that of increasing po as shown in Figure 7a-d, in terms of the improving chances of survival for the corresponding species. It should be noted that the double mushroom curves of Figure 8b and c are similar to those observed by Alhumaizi and Abasaeed33 for competition of parent and mutant autocatalators. The robustness of the mutant is inherited from the parent species, as continuous

3638


feeding of the latter provides a permanent avenue for production of the former. 4. Dynamics of Competition In the previous sections, the effect of the competition of an invading species with a robust host species on the steady-state and bifurcation structure of the system was investigated. This steady-state analysis showed that, whenever the invading species was present in the system, new steady states emerged that reflected greatly reduced levels of the host-species population. Whether the system migrates to one of these new steady states when subjected to an invasion or whether it successfully recovers its original state depends strongly on the nature of the invading species and the manner in which it is introduced. The dynamic behavior of this system is analyzed in this section to gain a better understanding of the response of the host-species population to the introduction of a second species. Specifically, we investigate the system’s response to a disturbance introduced as a rectangular concentration pulse of invading species in the feed. In the following case studies, the duration of the pulse in qo is fixed at 2 time units, and its magnitude is varied to test the robustness of the host species P. We also explore effective methods of combating these types of invasions to restore the original state of the system. An open-loop control strategy that adjusts the dimensionless residence time, T, of the reactor to a value outside the range of the stable steady states of the invading species Q is analyzed. This strategy can be implemented either by increasing T to starve to death or by decreasing it to wash out the invading species Q. In this section, we investigate the dynamic behavior of the system for the parameter sets of Figures 2b and 7d using MATLAB32 simulations. Recall that, in Figure 2b, the mushroom curve in the bifurcation diagram belongs to the single-species steady states, i.e., q ) 0, whereas the isola belongs to the twospecies steady states. An initial condition is selected on the upper steady state of the mushroom curve in p (lower steady state in r) for T ) 7.040 47, which is located in the multiple steady-state region of the mushroom curve and coincides with the range of stable steady states of the isola. If we perturb the system with a pulse qo ) 0.4 in the feed for t ∈ (10, 12), the system rejects the invasion on its own and stabilizes back at its initial steady state after a brief transient upset (Figure 9). In this case, the host species is robust enough to resist a perturbation of this magnitude. Increasing the size of the pulse of qo from 0.4 to 2.0, however, permanently disturbs the system as it perturbs the system from the upper stable steady state in p to the lower one on the mushroom curve in the (T, p) plane (Figure 10). Note that no Q is left in the system in both cases, suggesting that the perturbation size is not large enough for Q to gain a foothold in the system. In the second case, however, the disturbance is large enough to consume a considerable amount of resource R and to suppress the presence of P in the system, shifting the steady state of the system from a high concentration to a low concentration of P. If the size of the perturbation is further increased to qo ) 4.0, then Q survives in the system, as the system’s dynamics converges onto the stable steady state on the isola in Figure 2b (Figure 11). One way to remove Q from the reactor is to change T in either the increasing or decreasing direction, so that the operating point is shifted out of the stable region of

Figure 9. Effect of a pulse disturbance, qo ) 0.4, given in the feed for t ∈ (10, 12) on the system dynamics when started from a steady state at T ) 7.040 47 and r ) 0.250 584, p ) 0.651 755, and q ) 0.

Figure 10. With the same initial conditions as in Figure 9, a pulse, qo ) 2.0, is introduced for t ∈ (10, 12), which brings the system to the new steady state at higher concentration r and lower concentration p, with no Q left in the system.

Figure 11. Pulse, qo ) 4.0, is introduced to the system for t ∈ (50, 52). To remove Q from the system, T is increased from 7.040 47 to 20 at t ) 200. To restore the system to its initial conditions, T is decreased from 20 to 7.040 47 at t ) 400.

the isola. Figure 11 depicts the case where T is increased from 7.040 47 to 20 at t ) 200 to starve Q to death, for which the reactor goes to the steady state located at the


Figure 12. In contrast to Figure 11, T is decreased to 0.08 at t ) 200 to remove Q. Then, T is increased from 0.08 to 20 at t ) 210 and decreased from 20 to 7.040 47 at t ) 300 to restore the system to its initial conditions.

Figure 13. From Figure 4d, the initial condition T ) 47.9788, r ) 0.796 072, p ) 0.075 708, and q ) 0, which falls on the limit cycle, is selected. A pulse, qo ) 4.0, is introduced for t ∈ (1000, 1002). T is reduced to 2 at t ) 1200 and increased to 47.9788 at t ) 1400 to restore the system to its initial conditions.

single steady-state region of the mushroom curve at T ) 20. To bring the system back to its initial steady state, T is decreased back to 7.040 47 at t ) 400. We can also attain the original operating point by accelerating the flow through the system (Figure 12), which corresponds to a decrease in the residence time from 7.040 47 to 0.08 at t ) 200. However, in this case, after Q is removed from the reactor, we cannot simply increase T from 0.08 to 7.040 47 to restore the system to its initial condition, because the system will be trapped at the lower stable steady state in p on the lower-residence-time branch of the mushroom curve in the (T, p) plane. Thus, T has to be increased to 20, for example, which is beyond the multiple-steady-state region, to ignite P and then decreased back to 7.040 47. Note that species P is neither starved to death nor washed out with species Q because of the constant supply of po through the feed that ensures its robustness in the system. The bifurcation diagrams of Figures 2-8 show various regions of operation at which branches of stable fixed points become unstable at Hopf bifurcation points. Under these operating conditions, the possibility that sustained oscillatory solutions will emerge exists. The

Figure 14. Phase plane (r, p) shows two limit cycles of Figure 13.

Figure 15. Bifurcation diagram generated by increasing dq in Figure 4d from 0.0240 to 0.15 keeping kq ) 64dq to shift the isola to lower residence time shows the intersection of the unstable steady-state regions of the single-species and two-species curves.

dynamics of interaction between the robust host and the nonrobust invader will thus be impacted by the presence of oscillations. This kind of dynamic behavior is investigated next. In Figure 7d, the mushroom curve in the bifurcation diagram represents the single-species steady states, corresponding to a single solution (point c in Figure 1), and the isola represents the two-species steady states. An initial condition is selected on the limit cycle of the single-species steady states at T ) 47.9788 for which the second species has a stable steady state on the isola. Figure 13 shows the limit-cycle response starting at the selected initial condition with r ) 0.796 072 and p ) 0.075 708. When this state is perturbed by a pulse of qo ) 4.0 for t ∈ (1000, 1002), the limit cycle disappears as the system stabilizes at the stable steady state of the isola. To remove Q, we reduce T from 47.9788 to 2 at t ) 1200 and then increase it from 2 back to 47.9788 at t ) 1400 to restore the system to its initial condition. The simulation results show that Q is eliminated from the system using this strategy but that r and p are brought to a different limit cycle on the periodic branch at T ) 47.9788 with a decreased amplitude of oscillation. The phase plane of r and p in Figure 14 illustrates this phenomenon, as the two limit cycles differ in size and location. This

3640


5. Conclusions

Figure 16. From the bifurcation diagram in Figure 15, the initial condition T ) 39.2117, r ) 0.524 211, p ) 0.119 414, and q ) 0 on the limit cycle of the single-species system is selected. After a pulse, qo ) 5.0, introduced for t ∈ (500, 502), the system dynamics migrates to a two-species limit cycle. To remove q and restore the initial conditions, T is reduced from 39.2117 to 0.2 and then increased back to 39.2117 at t ) 1200 and 1400, respectively.

The system studied is an isothermal CSTR with two cubic autocatalytic reactions where both autocatalysts exhibit decay. The host species P is robust in the sense that it is constantly supplied to the reactor, whereas the invading species Q is introduced only as a pulse disturbance in the feed. The system has up to three steady states with q ) 0 and up to four additional steady states with q * 0. The steady states with q ) 0 form either an isola and an single steady-state curve or a mushroom on the bifurcation diagram, with T being the bifurcation parameter. Steady states with q * 0 form either one or two isolas, often with mushroom-like distortions. When the ratio kq/dq is kept constant, the major effect of changing the death rate of the invading species, dq, is on the location(s) of the isola(s). We observed that the robustness of the host species is decreased dramatically by the presence of the invading species when their death rates (dp and dq) are comparable. As the feed concentration of species P, po, is increased, the distortion of the isolas in the (T, r) and (T, p) planes is increased, and the isolas in the (T, p) plane grow bigger. When species Q is also fed continuously to the reactor, the isolas of the two-species steady states open up to give a mushroom. We have designed and tested an open-loop control scheme where the residence time, T, of the reactor is manipulated outside the stable steady-state region of the invading-species isola to eliminate this species from the reactor. This can be done either by increasing T to starve Q to death or by decreasing T to wash Q out of the system. The development of a closed-loop implementation of these strategies will be reported elsewhere.35 Literature Cited

Figure 17. Phase plane (r, p) shows two limit cycles of Figure 16.

behavior is a result of multiplicity on the periodic branch of solutions. Next, we investigate the system at an initial condition located in the unstable region of the mushroom that coincides with the unstable region of the isola as well. The bifurcation diagram of this case is shown in Figure 15. The simulation result in Figure 16 shows that the initial condition, selected at T ) 39.2117, where r ) 0.524 211, p ) 0.119 414, and q ) 0, leads to the limit cycle emanating from single-species steady-state branch. After the introduction of a pulse qo ) 5.0 for t ∈ (500, 502), the system is brought to a limit cycle of coexistence of the two species, where r, p, and q show oscillatory behavior with different amplitudes. The oscillation of p is hardly seen in Figure 16 because its amplitude is so small. To remove Q and restore the initial limit cycle, T is reduced from 39.2117 to 0.2 and then increased back to 39.2117 at t ) 1200 and 1400, respectively. The simulation result in Figure 16 shows that the initial limit cycle can be restored. The phase plane of r and p in Figure 17 illustrates the switch back and forth between the two limit cycles.

(1) Lin, K. F. Concentration Multiplicity and Stability for Autocatalytic Reactions in a CSTR. Can. J. Chem. Eng. 1979, 57, 476. (2) Lin, K. F. Multiplicity, Stability and Dynamics for Isothermal Autocatalytic Reactions in CSTR. Chem. Eng. Sci. 1981, 36, 1447. (3) Varma, A.; DeVera, A. L. Dynamics of Selectivity Reactions in Isothermal CSTRs. Chem. Eng. Sci. 1979, 34, 1377. (4) Heinrichs, M.; Schneider, F. W. On the Approach to Steady States of Reacting Systems in CSTR. Ber. Bunsen-Ges. Phys. Chem. 1980, 84, 857. (5) Escher, C.; Ross, J. Multiple Ranges of Flow Rate Bistability and Limit Cycles for Schlo¨gl’s Mechanism in a CSTR. J. Chem. Phys. 1983, 79, 3673. (6) Gray, P.; Scott, S. K. Absence of Restrictions on Permissible Paths for Changing the Efficiencies of Operation of Flow Reactors. J. Chem. Phys. 1983, 87, 1835. (7) Gray, P.; Scott, S. K. Autocatalytic Reactions in the Isothermal CSTR. Chem. Eng. Sci. 1983, 38, 29. (8) Gray, P.; Scott, S. K. Autocatalytic Reactions in the Isothermal CSTR. Chem. Eng. Sci. 1984, 39, 1087. (9) Scott, S. K. Reversible Autocatalytic Reactions in an Isothermal CSTR: Multiplicity, Stability and Relaxation Times. Chem. Eng. Sci. 1983, 38, 1701. (10) Li, R. S.; Li, H. Isolas, Mushrooms and Other Forms of Multistability in Isothermal Bimolecular Reacting Systems. Chem. Eng. Sci. 1989, 44, 2995. (11) Li, R. S. CSTR with Two Inflows of Reactants: A Versatile Tool for Study of Bifurcation in Chemical Systems. Chem. Eng. Sci. 1994, 49, 2029. (12) Focant, S.; Gallay, Th. Existence and stability of propagating fronts for an autocatalytic reaction-diffusion system. Physica D 1998, 120, 346.

Ind. Eng. Chem. Res., Vol. 41, No. 15, 2002 3641 (13) Stadler, B. M. R.; Stadler, F. P. Small Autocatalytic Reaction Networks III. Monotone Growth Functions. Bull. Math. Biol. 1991, 53, 469. (14) Stadler, F. P.; Schnabl, W.; Forst, C. V.; Schuster, P. Dynamics of Small Autocatalytic Reaction Networks II. Replication, Mutation and Catalysis. Bull. Math. Biol. 1995, 57, 21. (15) Murray, J. D. Mathematical Biology; Springer-Verlag: Berlin, 1993. (16) Varetto, L. Studying Artificial Life with a Molecular Automaton. J. Theor. Biol. 1998, 193, 257. (17) Kauffman, S. Nature 1996, 382, 496. (18) Orgel, L. E. Molecular replication. Nature 1992, 358, 203. (19) Kauffman, S. At Home in the Universe: The Search for Laws of Self-Organization and Complexity; Oxford University Press: New York, 1995. (20) Sharov, A. A. Signs and Values. In Proceedings of the IEEE ISIC/CIRA/ISAS Joint Conference; IEEE Press: Piscataway, NJ, 1998; p 677. (21) Wills, P. R.; Kauffman, S.; Stadler, B. M. R.; Stadler, P. F. Selection Dynamics in Autocatalytic Systems: Templates Replicating Through Binary Ligation. Bull. Math. Biol. 1998, 60, 1073. (22) da Silva, L.; Mundim, K. C.; Tsallis, C. Effects of crosslinks on the autocatalytic polymerization of RNA-like chains. Physica A 1998, 259, 415. (23) Schuster, P.; Fontana, W. Chance and necessity in evolution: Lessons from RNA. Physica D 1999, 133, 427. (24) Kauffman, S. Escaping the Red Queen Effect. McKinsey Q. 1995, 118. (25) Blomberg, C. On the Appearance of Function and Organisation in the Origin of Life. J. Theor. Biol. 1997, 187, 541. (26) Hofbauer, J.; Weibull, J. W. Evolutionary Selection against Dominated Strategies. J. Econ. Theory 1996, 71, 558.

(27) Hansen, S. R.; Hubbell, S. P. Single Nutrient Microbial Competition: Qualitative Agreement between Experimental and Theoretically Forecast Outcomes. Science 1980, 207, 1491. (28) Fredrickson, A. G.; Stephanopoulos, G. Microbial Competition. Science 1981, 213, 972. (29) 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. (30) Birol, I.; Teymour, F. Statics and Dynamics of Multiple Cubic Autocatalytic Reactions. Physica D 2000, 144, 279. (31) Kuznetsov, Yu. A.; Levitin, V. V. CONTENT: A Multiplatform Environment for Analyzing Dynamical Systems; Dynamical Systems Laboratory, Centrum voor Wiskunde en Informatica: Amsterdam, 1996. (32) MATLAB, version 5.3; The MathWorks, Inc.: Natick, MA, 1999. (33) Alhumaizi, K. I.; Abasaeed, A. E. On Mutating Autocatalytic Reactions in a CSTR. Chem. Eng. Sci. 2000, 55, 3919. (34) Birol, I.; Parulekar, S. J.; Teymour, F. Effect of Environment Partitioning on the Survival and Coexistence of Autocatalytic Replicators. Phys. Rev. E 2002, manuscript submitted. (35) Chaivorapoj, W.; Birol, I˙ .; C¸ inar, A.; Teymour, F. Feedback Control of CSTR with Competing Autocatalators. Ind. Eng. Chem. Res. 2002, manuscript submitted.

Received for review April 2, 2001 Revised manuscript received February 2, 2002 Accepted April 19, 2002 IE010302P

Competition between Robust and Nonrobust Autocatalytic Replicators

Recommend Documents