Analysis of a continuously stirred two tank reactor cascade with Haldane kinetics

Biological reactors are employed in industrial applications to break down organic waste. We view the cascade of two open loop continuously stirred tank reactors with Haldane growth function as chemostats with bacterial inputs. A function of some of the reactor parameters is derived, the sign of which determines the maximum number of critical points a reactor can have. This allows us to determine the parameter combinations which ensure a reactor has only a single critical point for all bacterial removal rates (dilution rate plus death rate). Where a simple condition on the above function is confirmed to hold, if the first reactor in a cascade only has a single critical point for all bacterial removal rates then, the next reactor will also only have a single critical point for all bacterial removal rates. A global stability result is also given for some of these cases. A simple proof is given for the local stability of critical points of a reactor with a general class of bacterial growth functions, bacteria and substrate input, and a death rate. For the special case where the first reactor has zero bacteria input, we compare a two reactor cascade with a single reactor under various conditions, long and short residence times, and different death rates. This follows the pattern of similar papers that considered cascades using the Monod and Contois growth functions.

called a Chemostat, have been studied extensively, cascades of reactors (reactor chambers in sequence) only gained prominence at a 1962 Prague Symposium with a paper entitled 'Multi-stage Continuous Culture' by Herbert (1964).
Herbert said "If one fermenter gives good results, two fermenters will give better results and three fermenters better still. This is sometimes true, but often false". He went on to note that for the purposes of producing quantities of cells (quantity of product, rather than reducing concentration of waste) it was seldom of practical benefit to use more than two reactors in a cascade.
The Chemostat however, being continuously stirred, is assumed to have the same concentrations of substrate and bacteria at every point in the reactor. Natural systems, though, have nutrient gradients and often more than one competing bacterial species. To model these properties, Lovitt and Wimpenny introduced a cascade of Chemostats, called a Gradostat, with two way flow between adjacent Chemostats, see Lovitt and Wimpenny (1979) and Lovitt and Wimpenny (1981). The two way flow between reactors discretely models the diffusion that would occur in a continuous medium. This approach has been generalised to consider arbitrary arrangements of connected Chemostats, see Rapaport (2018). A detailed treatment of the Chemostat and Gradostat can be found in the texts by Smith and Waltman (1995) and Harmand et al. (2017) and the extensive literature referenced therein. In this paper we restrict ourselves to a reactor chain, where the output of one reactor simply feeds into the next reactor.
While various reactor parameters affect the output of a reactor cascade, probably the most significant is the growth function, that function describing the reaction between bacteria and substrate throughout the reaction. A range of growth functions is well known in the literature and studies of cascades of some of these have been reported. Thus, Sidhu et al. (2015) have studied a reactor cascade which uses the Monod growth function, while Nelson and Holder (2009) have done so for a cascade using the Contois growth function. Dramé et al. (2006), looking at multiple steady states in interconnected biological systems, considered a cascade of reactors with a general growth function, which included the Haldane growth function as a special case.
The Haldane growth function was introduced to account for the observation that higher concentrations of substrate, past a certain point, can suppress the growth rate of bacteria (Haldane 1930). The Haldane model has been found to fit a wide range of reactor data for bacterial growth suppressed by toxic substrates (Muloiwaa et al. 2020). The wide applicability of this model makes its study important.
We will also include a death rate in our reactors. We note that Herbert (1958) has suggested that for long residence times the discrepancy between the predictions of the Monod model and the observed values can be corrected by the inclusion of a bacterial death rate. The inclusion of a death rate in the models considered here allows us to explore the potential impact of a death rate on a reactor. A death rate can also be included by allowing the removal rate to differ for each bacterial species, as is done in Li (1998).
The main difficulty one encounters when studying a cascade of reactors is a lack of simple analytic expressions for the substrate and bacteria concentrations at critical points, when a reactor has both substrate and bacterial inputs. For this reason we choose to study the simplest case of a cascade without recycle.
Finally, we note that optimization of the output of the cascade is not considered. This paper contributes to the literature by looking at a range of properties of a two reactor cascade with Haldane growth function. We view the cascade of two open loop continuously stirred tank reactors with Haldane growth function as chemostats with bacterial inputs. A function of some of the reactor parameters is derived, the sign of which determines the maximum number of critical points a reactors can have. This allows us to determine the parameter combinations which ensures a reactor has only a single critical point for all bacterial removal rates (dilution rate plus death rate).
Where a simple condition on the above function is confirmed to hold, if the first reactor in a cascade only has a single critical point for all bacterial removal rates then, the next reactor will also only have a single critical point for all bacterial removal rates. A global stability result is also given for some of these cases.
A simple proof is given for the local stability of critical points of a reactor with a general class of bacterial growth functions, bacteria and substrate input, and a death rate.
For the special case where the first reactor has zero bacteria input, we compare a two reactor cascade to a single reactor under various conditions, long and short residence times, and different death rates. This follows the pattern of similar papers that considered cascades using the Monod and Contois growth functions.
The structure of this paper is as follows. Section 2 presents a dimensional system of equations to model the cascade of two reactors with Haldane growth function and bacterial input. This system is then converted to a dimensionless autonomous system of equations.
Section 3 considers steady state solutions of this system by first presenting a brief review, Sect. 3.1, of the critical points and their stability for Reactor 1 with no bacterial input. Section 3.2 looks at the critical points, Sect. 3.2.1, and stability, Sect. 3.2.2, of Reactor 2. It is important to note that this stability is established using linear stability theory, but without the use of explicit analytic expressions for the critical points. We show that the stability of a critical point is related to the slope of r 2 , a function generalising the growth function, that we will define. This relationship is given in Theorem 1.
In Sect. 3.2.3, we introduce , a function of K , the single parameter used in the Haldane growth function, and the input substrate and bacteria concentrations and in Theorem 2 show that the sign of gives a sufficient condition for the reactor to have a single critical point for all bacteria removal rates, which in addition is locally asymptotically stable. In Sect. 3.2.4, we look at the number of critical points in a single reactor and relate this to the parameters of the reactor. In Sect. 3.2.5, using this function, we see that for certain combinations of reactor parameters in the first reactor, where it has a single critical point for all bacterial removal rates, if the = 0 curve for the given K is a non-decreasing function of bacteria concentration x in the domain of interest, it is possible to conclude that the second reactor would also only have a single critical point for all bacterial removal rates. In Sect. 3.2.6, we provide a global stability result for the first reactor in the cascade with some conditions on the inputs. In Sect. 3.2.7, we provide a global stability result for the second reactor and two reactor cascade.
Section 4.1 considers the dependence of washout for Reactor 1, with no bacterial input, on the input substrate concentration, death rate and hydraulic residence time. Section 4.2 compares the performance of a two reactor cascade (with no bacterial input to the first reactor) to that of a single reactor (with no bacterial input) for large residence times, but without a death rate. Section 4.3 then compares the performance of a two reactor cascade(with no bacterial input to the first reactor) to that of a single reactor (with no bacterial input) for residence times nearer to washout but with a few different death rates.
Section 5 presents a discussion of these results and shows how Theorem 1 can be used to extend some of the results given in Dramé et al. (2006).
In what follows the reader is asked be aware that wherever f (s) appears the statement or equation will always be true for the Haldane growth function. However we note that Eqs. (10) and (11), (19) to (24), and Theorem 1, will also be true for any f (s) which is differentiable, passes through the origin and is non-negative and bounded on s > 0. That is Theorem 1 Fig. 1 Two reactor cascade; S in , S 1 , and S 2 are substrate concentrations, X in , X 1 , and X 2 are bacterial concentrations, Q is the common volume flow rate, and V 1 and V 2 are the two tank volumes is also valid for functions f (s) with multiple peaks and any number of critical points and includes the Haldane growth function as a special case.

Model
The model considered here is a cascade of two continuously stirred tank reactors without recycle, shown in Fig. 1. The first reactor, Reactor 1, has a continuous inflow of substrate with concentration S in and bacteria with concentration X in , and an outflow of both substrate and bacteria concentrations, S 1 , X 1 respectively. The outflow of Reactor 1 is then fed directly into Reactor 2, with its subsequent outflow of substrate and bacteria concentrations, S 2 , X 2 respectively. The volume flow rate, Q, is assumed to be the same for both reactors to allow for continuous operation. However the volumes of the two reactor tanks, V 1 , V 2 may differ, giving different dilution rates, D 1 , D 2 . Note that the hydraulic residence times, the mean time a substrate or bacterium particle is in a reactor, are the inverses of the dilution rates, 1/D 1 , 1/D 2 .
As the bacterial species is assumed to be the same in both tanks, the bacterial death rates and growth rate functions, including the maximum specific growth rate, will be the same for both reactors.
Here, S 1 (T ), S 2 (T ) and X 1 (T ), X 2 (T ) are the substrate and bacterial concentrations respectively, at time T , where the subscripts 1 and 2 indicate Reactor 1 and Reactor 2 respectively. The concentrations of substrate and bacteria entering Reactor 1 are assumed to be constant.
The dilution rates for each reactor are given by D 1 = Q/V 1 , D 2 = Q/V 2 , where V 1 , V 2 are the volumes of the two reactor tanks.
F(S) is the common growth rate function with parameters, K s , the substrate concentration at half maximum specific growth rate and, K i , the inhibition constant. The maximum specific growth rate, and death rate for the bacteria are M and K d respectively. The common yield coefficient, Y , assumed constant, relates the bacterial growth rate to the substrate utilization rate. Here, we assume all parameters to be positive.

The dimensionless model
We will work with a dimensionless version of the problem, by introducing the dimensionless variables s n (t), x n (t), n = 1, 2 and t, and the dimensionless parameters, K , k d , s in , x in and d n , In terms of these dimensionless quantities, the equations that govern the two reactor cascade become, where We note that the dimensionless parameters that are involved with the removal of bacteria from a reactor are the dilution rates d n , n = 1, 2 and the death rate k d . The sum, b n = d n + k d , n = 1, 2, will become important when we find the critical points for each reactor. As noted above, the dimensionless hydraulic residence times, τ n = 1/d n , n = 1, 2, are the inverses of the dimensionless dilution rates.

Steady state solutions
A steady state, or Critical Point(CP), exists when all of the state variables, for either a single reactor, s 1 , x 1 , or the cascade s 1 , x 1 , s 2 , x 2 , no longer vary in time. In finding steady states and their stability we will use the theory for autonomous systems, where the single reactor is governed by (8), (9), (12), and the cascade by (8)-(12).
As we assume s in and x in to be constant, and the system is open loop (there is no feedback from either reactor), the equations governing Reactor 1 will be autonomous. Similarly if Reactor 1 is operating at a steady state, its output bacterial and substrate concentrations will be constant, hence the input bacterial and substrate concentrations for Reactor 2 will be constant, and the equations governing Reactor 2 will be autonomous.
The CPs for Reactor 1, (s * 1 , x * 1 ), are found by setting the right hand sides of (8) and (9) to zero and solving.
The CPs are not all stable however. A CP is Asymptotically Stable (AS), when, after a small displacement from the CP, the system returns to the CP-a local property.
We are only interested in critical points where all of the four concentrations are nonnegative, that is we require s * 1 , x * 1 , s * 2 , x * 2 ≥ 0. In addition, since substrate cannot be created in a reactor, we require s * 1 ≤ s in and s * 2 ≤ s * 1 . We note that Reactor 1, with zero bacterial input, has been extensively studied by many authors Lara-Cisneros et al. (2012), Smith and Waltman (1995), Ajbar andAlhumaizi (2011), Harmand et al. (2017), and many others, and so here we will only give a summary of results that are needed in subsequent sections. We note however that not all analyses include a death rate. Moreover, our analysis of Reactor 2 (and Reactor 1 with a bacterial input) is further complicated by an inability to find simple analytic expressions for the critical points. This in turn makes the application of linear stability analysis in these cases less straightforward.

Critical points and stability of Reactor 1 when x in = 0
The critical points of the first reactor (s * 1 , x * 1 ), when x in = 0, are solutions of (8) and (9) where the right hand sides are set to zero, that is, they are solutions of the system of two equations We first find x 1 in terms of s 1 by adding (13) and (14) and rearranging to give where Substituting for x 1 from Eq. (15) into Eq. (13) results in the equation for s 1 Note that s * 1 = s in , the washout solution, exists for every b 1 > 0. We see from Fig where the value of, b 1 > 0, is independent of f (s 1 ). These solutions can be seen graphically in Fig. 2 as the points of intersection of the functions b 1 and f (s 1 ).
At the CP (s * 1 , x * 1 ), we can interpret the condition b 1 = f (s * 1 ) physically as a balance between the growth rate of live bacteria, f (s * 1 ), and the rate at which live bacteria is removed from the reactor, b 1 = d 1 + k d .
Recalling that for a CP (s * 1 , x * 1 ), we require s * 1 ≤ s in , we see the number of CPs which satisfy this, where the s * 1 are solutions of b 1 = f (s 1 ), depends on both K and s in . When s in ≤ 1/ √ K there is only one CP which is also AS, while for s in > 1/ √ K , there are two CPs, the first of which is AS, while the second is unstable. Figure 2 displays the stability of all such points, where AS CPs are indicated by solid red lines and unstable CPs by dashed red lines. The stability of the CPs (s * 1 , x * 1 ) for s * 1 mentioned above are well known (Dochain and Vanrolleghem 2001) and easily calculated directly using linear stability theory. Note , there are three s * 1 values for each b 1 , one stable and one unstable being solutions of (18) and the stable solution s * 1 = s in . It is possible for Reactor 1 to jump from one stable state to the other if the system is disturbed.
In Fig. 2 b 1 is assumed to take on values in the interval (0, ∞), if however the death rate k d > 0, then b 1 would be restricted to the interval (k d , ∞).

Critical points and stability of Reactor 2
In Sects. 3.2.1 and 3.2.2 below we consider the critical points of Reactor 2 and their stability. We remind the reader however that for constant inputs, s * 1 and x * 1 , to Reactor 2, (10) and (11) can be treated as an independent autonomous system of differential equations. We note that the analysis that follows could be repeated for the case of Reactor 1 with x in > 0.

Critical points of Reactor 2
The critical points of Reactor 2, (s * 2 , x * 2 ), are solutions of (10) and (11) where the right hand sides are set to zero, s 1 = s * 1 and x 1 = x * 1 , that is, they are solutions of the system of two equations where s * 1 and x * 1 are constant. Using a procedure analogous to that used in Reactor 1 with x in = 0 we first find x 2 in terms of s 2 by adding (19) and (20) and rearranging to give where Substituting (21) into (19) and rearranging, the s * 2 values of the CPs, (s * 2 , x * 2 ), are then solutions of where Note the choice of the two functions, b 2 and r 2 (s 2 ), in (23). The constant function b 2 contains the information about the removal of bacteria from the reactor and can be considered to be a bacteria dependent removal rate of the form considered by Li (1998). While r 2 (s 2 ), considered to be a generalization of f (s 2 ), contains the information about the growth of the bacteria in the reactor, where bacterial growth depends on the properties of the bacteria through f (s 2 ) and the concentrations of the input bacteria x * 1 and substrate s * 1 . We will see later that the shape of r 2 (s 2 ) alone determines the maximum number of critical points that are possible for Reactor 2. This is a natural extension of the case for Reactor 1 with x in = 0 where the maximum possible number of critical points depended on K and s in . The number of critical points will be independent of b 2 for some combinations of x * 1 , s * 1 and K as we will see. Equation (23) can also be rearranged as the cubic equation in s 2 , The CPs for Reactor 2 are given by obtaining the real positive solutions of (23), for s * 2 ∈ (0, s * 1 ), then using (21) to find the corresponding x * 2 values. From (23) we note that the critical points are the points of intersection of the constant function b 2 and r 2 (s 2 ). In addition since the value of b 2 is independent of s * 1 , x * 1 and K , it can be varied independently of r 2 (s 2 ). Changing the flow rate Q, which can be done by an operator, or the reactor tank volume V 2 , which can be incorporated into the design, will change the value of b 2 .
The requirement that b 2 = r 2 (s 2 ) at a CP can be interpreted physically as the balance of the rate of removal of live bacteria from Reactor 2, given by b 2 , with the rate of replenishment of live bacteria, given by r 2 (s 2 ). Typical plots of r 2 (s 2 ) are shown in Fig. 3. To see that Reactor 2 will never washout, note that r 2 (s 2 ) has a vertical asymptote at s 2 = s * 1 , so s * 2 → s * 1 as b 2 → ∞, but never reaches it. Although Reactor 2 never washes out, if s * 2 is close to s * 1 , very little substrate will be converted.

Stability of Reactor 2
One of the contributions of this paper is to present a simple proof of the stability of the CPs of Reactor 2, which include both a death rate and influent bacteria, using linear stability analysis in the absence of analytic expressions for the solutions of (23).
The stability of the CPs is determined by looking at the eigenvalues of the linearised system of equations evaluated at each of the critical points (s * 2 , x * 2 ). The Jacobian matrix of the right hand sides of (10)-(11), at the critical point is while the determinant and trace are given by, The eigenvalues are then given by When tr J (s * 2 , x * 2 ) < 0 and det J (s * 2 , x * 2 ) > 0 both eigenvalues will have negative real parts, so the CP will be AS. The critical point, (s * 2 , x * 2 ), will be a saddle point if the eigenvalues are real and of opposite sign.

Theorem 1 For all CPs
. Three cases exist: is a saddle node bifurcation point. All proofs can be found in Appendix.
The determinant has the same sign as r 2 (s * 2 ) and is zero at points s * 2 where r 2 (s * 2 ) = 0, being points where the stability of the critical points changes.
We remind the reader that this result, although provided in the context of a second reactor in a cascade, applies equally well to a single reactor with arbitrary but constant feed concentrations of substrate, s * 1 , and bacteria, x * 1 . This shape of r 2 (s 2 ), and hence the number of possible CPs for a given b 2 , is determined by the values of the input concentrations, x * 1 , s * 1 , and K alone. This echoes the case for Reactor 1 with x in = 0, where the number of CPs depended only on the parameters s in and K .
In summary, we have used local stability analysis to determine the stability of all the CPs, where r 2 (s * 2 ) > 0, are stable nodes, points where r 2 (s * 2 ) < 0, are saddle points, and points where r 2 (s * 2 ) = 0, are saddle-node bifurcation points.

Determining the nature of r 2 (s 2 ) for Haldane
The existence of bifurcations in the second reactor can be established using singularity theory (Ajbar and Alhumaizi 2011), which gives a set of six equations that the two saddle point bifurcation points must satisfy. Because these conditions refer only to the bifurcation points themselves they are functions of s * 1 , x * 1 , K and b 2 . However it is clear from Fig. 3 that the existence of bifurcations depends only on the shape of r 2 (s 2 ) and hence only on the values of s * 1 , x * 1 and K . The Theorem below gives a single function, , of s * 1 > 0, x * 1 > 0 and K > 0, the sign of which can establish the existence of bifurcation points. In essence this condition will tell us if the function r 2 (s 2 ) is monotonically increasing or not.

Theorem 2
The function r 2 (s 2 ), on the interval (0, s * 1 ), has: no real roots if < 0; a real double root if = 0; and two distinct real roots if > 0, where 1 , x * 1 > 0, and r 2 (s 2 ) given by (24). Figure 4 plots the = 0 boundary lines, in the positive quadrant, s * 1 ≥ 0, x * 1 ≥ 0, for a few different K values. For each line, > 0 above the line and < 0 below the line. We note that each curve cuts the s * 1 axis at s * 1 = 1/ √ K . Note that a plot over a limited domain is not sufficient to establish behaviour over a larger domain, however if a reactor has input substrate and bacterial concentrations of s * 1 and x * 1 respectively then it is sufficient to plot = 0 over the region (s 2 , x 2 ) ∈ (0, s * 1 )×(0, s * 1 + x * 1 ), for the given K . In what follows, and particularly in Sect. 3.2.5 below, we assume that the = 0 boundary, viewed as, s 1 a function of x 1 , examples of which are given in Fig. 4, for the given K , has been confirmed to be a non-decreasing function of x 1 over the domain of interest. Plots of these functions over extended ranges for the three variables have all been shown to be monotonically increasing.

Number of critical points in terms of system parameters
The number of critical points for a single reactor with input of both substrate and bacteria is considered here. Although the discussion is for Reactor 2, the argument is true for any reactor with Haldane growth function, and inputs of both substrate and bacteria.
We note that ≤ 0 is a sufficient condition for b 2 = r 2 (s 2 ) to have a single solution on (0, s * 1 ) for every b 2 ≥ k d , and hence for the reactor to have a single CP for all b 2 ≥ k d . This is illustrated by looking at the blue curve in Fig. 3.
If however > 0, then depending on the value of b 2 , Reactor 2 can have either a single CP, or three CPs. Where > 0, b 2 = r 2 (s 2 ) has three solutions (the reactor has three CPs) only for b 2 ∈ [b 2 min , b 2 max ], where b 2 max is the value of the local maximum of r 2 (s 2 ) on (0, s * 1 ) and b 2 min is the value of the local minimum of r 2 (s 2 ) on (0, s * 1 ). There are no easy analytic expressions for b 2 min and b 2 max unfortunately, however it is straightforward to numerically find the roots of r 2 (s 2 ) and hence the interval [b 2 min , b 2 max ].
Where b 2 has a value outside of the interval [b 2 min , b 2 max ], b 2 = r 2 (s 2 ) will have a single solution and the reactor a single CP. We note that where there are three CPs, if the system is disturbed the system may jump from one stable state to the other.

Single critical points in both reactors of a cascade
For a cascade of two reactors, where the first reactor has s in > 0, x in > 0 and K such that < 0 then the point (s in , x in ) will lie below the = 0 line for that K . The reactor will then have only a single CP for any bacterial removal rate.
Note that we do not require s in < 1/ √ K to ensure < 0. Also the bacterial removal rate can, in principle, be given any positive value greater than the death rate (by changing the tank volume), so the single critical point can, in principle, depending on the bacterial removal rate, have an s * 1 value at any point in the interval (0, s in ). This means it is possible for some choices of s in and bacteria removal rate that f (s * 1 ) < 0 when < 0. The input to the second reactor, (s * 1 , x * 1 ), will have s * 1 < s in and x * 1 > x in , so the input concentrations to the second reactor, (s * 1 , x * 1 ), will also lie below the = 0 line for that K . The second reactor will then have only a single CP for any bacterial removal rate.
We note that being independent of the bacterial removal rate means the result does not depend on the reactor tank volumes. This argument could be repeated for a third reactor and indeed any number of subsequent reactors.
For a cascade of two reactors where the first reactor has no bacterial input, x in = 0, and is not washed out, its critical point will have s * 1 < 1/ √ K and x * 1 > 0. Following a similar argument to that above, the second reactor will have < 0, and will have a single CP for any bacterial removal rate. The argument above tells us reactors past the second will all have a single CP for all bacterial removal rates.

Global stability of Reactor 1 with s in < 1/ √ K and x in ≥ 0
Here we consider the global stability of a single reactor, being Reactor 1, with Haldane growth function, and inputs of substrate and bacteria having constant concentrations, s in < 1/ √ K and x in ≥ 0 respectively. Where x in = 0, we also require x 1 (0) > 0, to ensure the reactor is not initially washed out.
In both cases, x in > 0 and x in = 0, a proof a global stability requires that only one AS CP exist and that the system evolve to this CP for all initial conditions. In the x in > 0 case, the conditions x in > 0 and s in < 1/ √ K are sufficient to ensure there is a single CP, which is AS, for any b 1 > k d ≥ 0. In the x in = 0 case a non-washout AS CP exists if b 1 < f (s in ) and x 1 (0) > 0.
With a slight abuse of notation (although the AS CPs for the two cases will be different in general) we will denote the AS CP for them both (s * 1 , x * 1 ). We note that neither case allows the possibility a cyclically chained orbit, so to establish global stability it remains to show that no limit cycle exists in the positive quadrant of s 1 > 0, x 1 > 0. The argument that follows will apply equally to both cases listed above.
As Region 2 does not contain a critical point it cannot contain a limit cycle (Jordan and Smith 1987 page 81).
To show that there is no limit cycle in Region 1 we will use Bendixson's Negative Criteria which says if the divergence of the system of the differential equations has a single sign in a connected domain then it does not contain a limit cycle (Jordan and Smith 1987 page 82). Using the change of variable ξ 1 = ln x 1 in the right hand sides of (8) and (9), and taking the divergence we get which is always negative in Region 1, where we note that ξ 1 can take any real value, so there are no limit cycles wholly contained in Region 1.
Finally we need to exclude the possibility that a limit cycle might exist in the combined Regions 1 and 2. To eliminate this possibility we note that if a limit cycle was partly in both regions ds 1 /dt would need to have a positive sign as it flowed from Region 1 to Region 2 across s 1 = s in , and a negative sign when it flowed back, but ds 1 /dt ≤ 0 for all x 1 on s 1 = s in . We conclude that there are no limit cycles in the positive quadrant and hence the single attracting critical point is a global point of attraction. We remind the reader that these arguments apply equally to the two cases, x in > 0, and x in = 0 with x 1 (0) > 0.
We note that if Reactor 1 had already reached a steady state, and the inputs to Reactor 2 were therefore constant, a similar argument could be used to establish the global stability of Reactor 2. However as the output concentrations of Reactor 1 vary in time for t > 0, the input concentrations to Reactor 2 will also vary in time, that is the system of equations governing Reactor 2 will be non-autonomous.
A global stability proof of the two reactor cascade therefore would, in addition, require a proof that the non-autonomous system governing Reactor 2 also approached a single AS CP for all initial conditions. As the input concentrations of substrate and bacteria for Reactor 2 approach the steady state concentrations for Reactor 1, the theory of asymptotically autonomous differential equations will apply. We consider this in the next section.

Global stability of Reactor 2
To show the two reactor cascade is globally asymptotically stable we need to show that both Reactor 1 and Reactor 2 each have a single attracting CP to which all solutions are attracted. We have already shown that, under the specified input concentrations, Reactor 1 is globally asymptotically stable. We now want to show that all solutions of the non-autonomous system governing Reactor 2 converge to a single attracting point.
Reactor 2 is governed by where s 1 (t) → s * 1 and x 1 (t) → x * 1 as t → ∞. We first note that the limiting case of (32) and (33) is where the system (32), (33) is said to be asymptotically autonomous with the limiting system (34), (35) if the right hand sides of (32), (33) tend to the right hand sides of (34), (35) as t → ∞, locally uniformly on s 2 ≥ 0, x 2 ≥ 0. As the right hand sides of (32), (33), (34) and (35) are all C 1 they are continuous and locally Lipschitz (Hirsch and Smale 1974 page 163), which is sufficient to conclude that the convergence is locally uniform. In addition we will assume that all solutions exist for all forward times.
Theorem 1.1 in Thieme (1994) (being a result from a paper by Marcus (1956)), tells us that the ω-limit set ω of a forward bounded solution (s(t), x(t)) to the non-autonomous system (32), (33) is non-empty, compact and connected. Moreover ω attracts (s(t), x(t)), that is, the distance between (s(t), x(t)) and ω tends to zero as t → ∞. Finally ω is invariant under the autonomous system (34), (35). In particular any point in ω lies on a full orbit of (34), (35), that is contained in ω.
We first note that (34), (35) is an autonomous system, and since s * 1 < 1/ √ K and x * 1 > 0, it has a single CP, (s * 2 , x * 2 ), which is AS. Also, by a similar argument to that given for Reactor 1, it contains no limit cycles. That is, (s * 2 , x * 2 ) is the only invariant point of this system. If we consider an arbitrary forward bounded solution of (32), (33), we can conclude, using Theorem 1.1 from (Thieme 1994), that it is attracted to its ω-limit set, which is non-empty, compact and connected. Finally ω is invariant under (34), (35). However the only invariant set of (34), (35) is the single AS CP, so the ω set must only contain the single AS CP (s * 2 , x * 2 ). Hence all forward bounded solutions to (32), (33) will be attracted to (s * 2 , x * 2 ). To prove the global stability of Reactor 2, and hence the cascade, we still need to show that every solution of (32), (33) is forward bounded. We start by showing that solutions of (34), (35), where k d = 0, are forward bounded.
Setting k d = 0 in (35) and adding (34) gives a differential equation for s 2 (t) + x 2 (t), whose solution is where are constants, and s 2 (0), x 2 (0) are initial concentrations. It is straightforward to see that A will bound solutions when B ≤ 0, while A + B will bound solutions when B > 0. So all solutions will be forward bounded when k d = 0.
To see that this remains true when k d > 0, consider the phase plane diagram for the k d = 0 case, with s * 1 and x * 1 fixed. The tangent vector at any point (s 2 , x 2 ) has components given by the right hand sides of (34), (35). We note that the s 2 component of the vector field is unchanged by the inclusion of k d > 0, while the x 2 component is decreased, that is the tangent vector will be inclined further towards the bounding s 2 axis than it does in the k d = 0 case. Since the solutions are determined by the vector field, adding a bacterial death rate will not lead to solutions going to infinity, so all the solutions of (34), (35) must also be forward bounded.
We now consider the forward boundedness of solutions of solutions of (32), (33). If a solution is to fail to be forward bounded it must escape to infinity either in finite time or infinite time. We first consider escape in finite time.
The argument above allows us to conclude that if s 1 (t) and x 1 (t) in (32), (33) were to be replaced by finite positive constants, the solution from the resulting system would be forward bounded. Since Reactor 1 is globally asymptotically stable, s 1 (t) and x 1 (t) never become infinite, so over a small finite time interval t, the solutions s 2 (t) and x 2 (t) to (32), (33) would be bounded, even if the values of s 1 (t) and x 1 (t) varied over the small time interval. This argument can be repeated for a finite number of time intervals t, so solutions to (32), (33) will be bounded for any finite time t.
To show that the solutions will not escape in infinite time we note that, (s 1 (t), x 1 (t)) can be made to approach (s * 1 , x * 1 ) as closely as we like by choosing a large enough finite time t. However since the convergence of the right hand sides of (32), (33) to (34), (35) is locally uniform the phase plane diagram for (32), (33) can be made locally as close as we like to that of (34), (35). As the solution paths are determined by the vector field depicted in the phase plane, solutions to (32), (33) could only escape to infinity if the solutions to (34), (35) were to escape to infinity. Since we know that the solutions to (34), (35) are forward bounded as t → ∞, the solutions to (32), (33) must also be forward bounded as t → ∞.
In conclusion, Reactor 2 will be globally asymptotically stable and hence the two reactor cascade will be globally asymptotically stable. We note in addition that, as the existence of the single AS CP for Reactor 2 holds for any b 2 > k d , the result will be independent of the relative volumes of the two reactors.

Two reactor cascade properties with x in = 0
In this section we compare the behaviour of a two reactor cascade, shown in Fig. 1, with that of a single reactor, where both the single reactor and the two reactor cascade have x in = 0. These results can then be compared with those of Sidhu et al. (2015) and Nelson and Holder (2009) for the Monod and Contois growth functions respectively.

Washout
We first consider washout in Reactor 1. There are two cases that exist, s in ≤ 1/ √ K , and s in > 1/ √ K . In both cases washout can occur if b 1 ≥ f (s in ). In the second case however it is also possible for Reactor 1 to be in a non-washout AS state when b 1 ∈ ( f (s in ), f (1/ √ K )). For b 1 in this range it is possible for the reactor to jump from a non-washout state to a washout state. Washout can be avoided in both cases by setting b 1 < f (s in ).
Note that for b 1 < f (s in ), the output substrate concentration, s * 1 , is the same for both cases, s in ≤ 1/ √ K , and s in > 1/ √ K . The percentage reduction in substrate concentration is therefore greater if s in > 1/ √ K . Although s * 1 is the same in both cases, the bacterial concentration, x * 1 , for the s in > 1/ √ K case is much higher. The condition, b 1 < f (s in ), can be rearranged to give a condition on the dimensionless hydraulic residence time, τ 1 , needed to avoid washout, We note τ 1 being positive requires that f (s in ) > k d . Figure 5 plots this minimum dimensionless hydraulic residence time needed to ensure washout will not occur, as a function of the input concentration, s in , for a number of K values, with k d = 0.01. Increasing the death rate will increase these minimum residence times. For a choice of K , residence times above the curve will ensure washout will not occur, while for times below the curve washout may occur.
For the cascade of two reactors, first note that if Reactor 1 washes out then Reactor 2 will not be fed any bacteria and can washout. If the volume of both reactors is the same then Reactor 1 washing out will cause Reactor 2 to washout. If V 2 > V 1 , then it is possible, given the increased residence time in Reactor 2, that it may not wash out even though Reactor 1 has.
If Reactor 1 is not washed out then Reactor 2 will have a bacterial input and not be washed out. In addition, as was shown above, Reactor 2 will have < 0 and have a a single CP.

Large residence times
In this section we will look at the output of the cascade for large residence times but restrict our consideration to the case of zero death rate.
To find asymptotic expansions of s * n and x * n , n = 1, 2, for large residence times (small d n values) for each reactor, one typically needs analytic expressions for s * n and x * n as functions of d n . Although this is possible for Reactor 1, there are no such simple expressions in the case of Reactor 2.
For Reactor 1, Eq. (18) gives s * 1 as a function of d 1 and substituting this into (15) gives x * 1 as a function of d 1 . Setting k d = 0 and finding the Taylor series expansions of s * 1 (d 1 ) and x * 1 (d 1 ) for small d 1 about d 1 = 0 gives: For Reactor 2 we do not have convenient expressions for s * 2 (d 2 ), x * 2 (d 2 ), however it is possible to find the Taylor series expansion for implicitly defined functions s * 2 (d 2 ) and x * 2 (d 2 ), when k d = 0, for small d 2 , about d 2 = 0, with inputs (s * 1 , x * 1 ), given by Although it is possible, in principle, to do the calculation for s * 2 with k d > 0, the expressions that result cover multiple pages of computer algebra which have no obvious simplification and hence yield no real insight.
To compare a single reactor with the cascade of two reactors we make the following assumptions. The two reactors have a common flow rate Q. The first reactor in the cascade has volume V 1 , dilution rate d 1 = Q/V 1 and residence time τ 1 = V 1 /Q. The second reactor has volume V 2 = RV 1 , dilution rate d 2 = Q/(RV 1 ) and residence time τ 2 = RV 1 /Q, where R = V 2 /V 1 > 0, is the ratio of the two tank volumes.
To make a fair comparison we choose the total residency time to be the same in both cases, hence the residency time for the single reactor will be τ = τ 1 + τ 2 = (R + 1)τ 1 . This corresponds to a single tank volume of (R + 1)V 1 with a dilution rate d 1 /(1 + R).
For a single reactor with residence time τ we set the dilution rate in (38) and (39) to To get the output of the cascade we first find the output of Reactor 1, using residence time τ 1 by setting the dilution rate in (38) and (39) to d 1 which gives On substituting (44) and (45) into (40) to (41), setting s * We see that for a single reactor, the output substrate concentration for very long residence times decreases as 1/((1 + R)τ 1 ) = 1/(τ 1 + τ 2 ), the inverse of the residence time, whereas for the two reactor cascade it decreases as 1/(τ 2 1 s in R) = 1/(τ 1 τ 2 s in ), that is, as the inverse of the product of the two residence times for the two reactors. The cascade of two reactors provides a clear increase in performance compared to a single reactor with the same total residence time. This is beneficial where reduction in the substrate concentration or an increase in product concentration is desired.

Output as a function of residence time
In this section we compare the output substrate concentrations of a two reactor cascade with a single reactor as a function of dimensionless hydraulic residency time for four bacterial death rates. As it is not possible to provide general formulae for the output substrate concentrations, we present some numerical examples. Here we assume the two reactors in the cascade have the same volume, V 1 = V 2 = V , with a common flow rate Q, and that the total residency times for the single reactor and the cascade are the same. However we no longer need to assume the death rate is zero. Figure 6 looks at the case where s in = 1/ √ K , and shows the ratio s out /s in for both a single reactor and a cascade of two reactor as a function of residence time for four bacterial death rates. √ K , and shows the ratio s out /s in for both a single reactor and a cascade of two reactors as a function of residence time for four bacterial death rates. Note that for some residence times, in both the single reactor and the cascade, there are two possible values for s out /s in values. This occurs because for these residence times the first reactor can be in one of two possible stable states, a washout state and a non-washout state. Although these have both been plotted Reactor 1 can only be in one of these two states at any one time.

Discussion
Here we have considered a cascade of two open loop CSTRs with Haldane growth rate function with a possible input of bacteria. Although many authors have considered a single reactor with a Haldane growth function, most do not include a bacterial input, and many omit a death rate. The model for reactors considered here, unless specified, includes both.
The introduction of a steady bacterial input to any reactor means the bacterial concentration can never go to zero, that is the reactor cannot washout. A bacterial input could be introduced either as a separate inflow to the tank as might happen when bacteria sediment from the outflow is recirculated or when the outflow from one tank, which is no longer sterile, containing substrate and bacteria, is fed into the next tank in a cascade.
For Reactor 2, r 2 (s 2 ), may be viewed as a generalisation of the growth rate function, f (s 2 ). We have shown that s * 2 values of the CPs, (s * 2 , x * 2 ), for this reactor occur at the intersections of b 2 = d 2 + k d and r 2 (s 2 ). Theorem 1 tells us the stability of a CP is given by the sign of the slope of r 2 (s 2 ) at such points of intersection. Importantly this result does not require the existence of explicit analytic expressions for s * 2 and x * 2 . Note that although the discussion is for Reactor 2, a similar argument would apply to any reactor with Haldane growth function and an input of substrate and bacteria.
As the proof of Theorem 1 does not use the functional form of the Haldane growth rate function, the proof would be valid for any growth rate function that was differentiable, passed through the origin and, was positive and bounded on s > 0. These last two conditions derive from biological considerations and require zero bacterial growth rate when there is no substrate, and a bounded positive growth rate when substrate is present.
We note that in our formulation the death rate is combined with the dilution rate in the function b n . Herbert (1958) commented that in the Monod case a death rate would be needed for large residence times to accurately predict the output concentrations of a reactor. Here we can see that for large residence times, when d n is small, k d can make up a significant fraction of b n = d n + k d , n = 1, 2, altering the location of the CP if ignored. It is reasonable to expect that Herbert's comment would be equally applicable to a larger class of growth rate functions.
For a reactor with Haldane growth function, Theorem 2 provides a function, , which when negative, tells us that the reactor will have a single CP for all bacteria removal rates. It follows that for a cascade of N reactors with Haldane growth function, (where the = 0 boundary, considered as substrate concentration a function of bacterial concentration, for that K has been shown to be a non-decreasing function on the domain of interest) if the first reactor has < 0, then none of the reactors will be washed out and each will have a single CP which will be AS. In addition, this will be independent of the volumes of the reactors. We note that if the input substrate concentration for a reactor in the cascade falls below 1/ √ K the single CP for that reactor and for reactors after that will possess global stability.
For an N reactor cascade, with Haldane growth function, where the first reactor has inputs of both substrate and bacteria, such that > 0, (where the = 0 boundary, considered as substrate concentration a function of bacterial concentration, for that K has been shown to be a non-decreasing function on the domain of interest) then it may have one or two AS CPs depending on the bacteria removal rate. However it will be the case, that for any reactor in the cascade, the output substrate concentration will be less than the input substrate concentration. This will result in a monotonic decrease in substrate concentration with each successive reactor in the cascade. If the cascade is long enough the substrate concentration will eventually fall below the = 0 curve (for the corresponding K value) and the CPs for each of the reactors from that point on will only have a single CP. A point also made by Dramé et al. (2006) in a similar context.
Section 4.3 clearly shows the accelerated reduction in substrate concentration obtained by a two reactor cascade compared to that of a single reactor. This parallels similar findings found for the Monod and Contois growth rate functions by (Sidhu et al. 2015), (Nelson and Holder 2009) respectively. Also illustrated in Figs. 6 and 7 is the decrease in this reduction of substrate with increasing death rate. This is to be expected given that increases in the death rate results in a reduction of bacteria available to convert substrate. We also note that the residence time needed to avoid washout increases with increasing death rate.
The work presented in this paper also generalises some of the work given in Dramé et al. (2006), by considering a general growth function F(S), which need only be once differentiable, and including a death rate. To see this consider an N reactor cascade where the first reactor has both a substrate and bacterial input, and the output of each reactor feeds directly into the next reactor. The cascade is governed by the system where i = 1, . . . , N and X 0 = X in and S 0 = S in are the concentrations of bacteria and substrate inputs into the first reactor. Now letting x i = X i /Y , d i = Q/V i , s i = S i , T = t, and F(S i ) = f (s i ), we get the system where i = 1, . . . , N , x 0 = x in and s 0 = s in . The critical points for each reactor can then be found recursively by solving the equations, where i = 1, . . . , N and s * i−1 and x * i−1 are constant inputs to Reactor i, (s * i−1 , x * i−1 ) being one of the stable operating points of the previous reactor, or the constant input concentrations in the case of the first reactor. Now setting b i = d i + k d , and adding (52) and (53), x i is given by = −K s 4 2 + 2K (s 1 + x 1 )s 3 2 + (−K s 1 (s 1 + x 1 ) + x 1 + 1)s 2 2 − 2s 1 s 2 + s 1 (s 1 + x 1 ) (s 2 − s 1 ) 2 (K s 2 2 + s 2 + 1) 2 .
Note that the coefficient of s 2 2 might be positive, negative, or zero, but in all cases the coefficients of the polynomial in s 2 will always only have one sign change. So applying Descartes' rule of signs to (60) we see that h(s 2 ) has exactly one real root in the interval (−∞, 0).
If s 1 ≥ x 1 , then the coefficient of s 3 2 will be negative or zero, so it does not matter what value the coefficient of s 2 2 has there will only be one sign change in the coefficients. Similarly if s 1 < x 1 the coefficient of s 3 2 will be positive, but s 1 < 5x 1 , so the coefficient of s 2 2 is also positive, again we only have one sign change. By Descartes' rule of signs, h(s 2 + s 1 ) has exactly one real root in the interval (0, ∞), so h(s 2 ) has exactly one real root in the interval (s 1 , ∞).
In addition h(0) = s 1 (s 1 + x 1 ) = 0 and h(s 1 ) = s 1 x 1 (K s 2 1 + s 1 + 1) = 0, so there are no real roots at the points 0 or s 1 . Since h(s 2 ) has at most four real roots there must be either no real roots or two real roots in the interval (0, s 1 ).
To decide on the number of distinct real roots in the interval (0, s 1 ) we look at the discriminant, 1 (Irving 2004), of h(s 2 ). If 1 > 0 then h(s 2 ) either has four distinct real roots or no real roots. As we have shown that h(s 2 ) always has at least two real roots then we must conclude that there are two distinct real roots in (0, s 1 ) if 1 > 0. We note that multiple roots can only occur when 1 = 0.
To sum up, there are no real roots of r 2 (s 2 ) in the interval (0, s 1 ) if < 0, a real double root if = 0, and two distinct real roots if > 0.