Model-based scale-up methodology for aerobic fed-batch bioprocesses: application to polyhydroxybutyrate (PHB) production

This work presents a general model-based methodology to scale-up fed-batch bioprocesses. The idea behind this approach is to establish a dynamics hierarchy, based on a model of the process, that allows the designer to determine the proper scale factors as well as at which point of the fed-batch the process should be scaled up. Here, concepts and tools of linear control theory, such as the singular value decomposition of the Hankel matrix, are exploited in the context of process design. The proposed scale-up methodology is first described in a bioprocesses general framework highlighting its main features, key variables and parameters. Then, it is applied to a polyhydroxybutyrate (PHB) fed-batch bioreactor and compared with three empirical criteria, that are traditionally employed to determine the scale factors of these processes, showing the usefulness and distinctive features of this proposal. Moreover, this methodology provides theoretical support to a frequently used empirical rule: scale-up aerobic bioreactors at constant volumetric oxygen transfer coefficient. Finally, similar process dynamic behavior and PHB production set at the laboratory scale are predicted at the new operating scale, while it is also determined that is rarely possible to reproduce similar dynamic behavior of the bioreactor using empirical scale-up criteria.


Latin symbols A x
Cross-sectional flow area m 2 ð Þ a 1 ; a 2 ; a 3 Empirical coefficients for P g calculation dimensionless ð Þ c 1 ; c 2 ; c 3 Empirical coefficients for X m Maximum residual cell concentration g L À Á X t Total biomass concentration (P þ X)

Introduction
The scale-up of fermentation reactors as a key biotechnological problem was first introduced in the 1940s when the industrial penicillin production began [41,43].Developments within this field continue until today looking for more efficient routes to reproduce at an industrial scale the same performance targets set at the laboratory.However, efforts have failed to establish a straightforward and general protocol to tackle this issue [4,17,26].Thus, a particular strategy is commonly implemented for each individual product, process and facility [41], requiring detailed understanding of the process, experience, intuition and good luck to successfully scale it up [17,39].
Although a plethora of works have emerged concerning the scale-up of biotechnological processes [10,18,20,26,35,41,43,47], they can be classified into three basic approaches: (1) experimental, (2) physical, and (3) fundamental [23,36,38,46].The first one is based on the designer extensive practice in scaling up processes [5,23,26,38,55], the second one uses dimensionless numbers, variables or relations to bond the same process at different scales [43,50] and the third one involves proper modeling of the process under consideration [17,29,30,38].Here, both experimental and physical approaches are known as traditional scale-up methods, which have been applied to a wide range of biotechnological processes since the 1960s [10,18,21,24,37,46,54].Nevertheless, latent-based multivariable scale-up methods (part of the experimental approach) have recently emerged as parallel solution to bioprocesses scale-up problem [14].This approach involves the use of available process data and fundamental knowledge from certain plant (e.g.pilot-scale unit) to obtain black-box models for a second plant (e.g.industrialscale unit) [49].This method pillar is the definition of the process latent variables that consider the relation of common variables between the plants as well as the relation between all the variables within each plant [13].Finally, the fundamental approach has only been explored in the last few years due to past limitations in solving complex models [4,17,26,38].
Bioprocesses scale-up has been commonly faced within the experimental approach, using empirical rules that intend to keep constant some key process parameter as the scale is increased [17,21,37].Here, the most popular parameters susceptible to be fixed are volumetric power consumption [26,54], volumetric mass transfer coefficient [10,47,52], impeller tip speed [26,54], Reynolds number [17,26], circulation time [55], mixing time [26,46] and dissolved oxygen concentration [17].Even though these criteria represent a simple scale-up procedure, each one is suitable for a narrow range of processes and operating conditions.Moreover, there is no agreement of which of them should be used for each individual problem [17,47], evidencing that the scale-up remains as an open question in the development of the biotechnology industry.
Therefore, seeking for a general methodology to scaleup fed-batch bioprocesses, this work takes advantage of the singular value decomposition (SVD) of the process Hankel matrix to determine an Output Impactability Index.This index enables both the quantification of each dynamics significance within the process and the establishment of the critical point of the fed-batch operating trajectory where the process should be scaled up.The methodology herein proposed is inspired on the procedure developed in [30] for chemical batch reactors, focusing on the particular characteristics of aerobic fed-batch bioprocesses and providing a theoretical justification, framed under the fundamental scale-up approach, for maintaining the same volumetric mass transfer coefficient throughout aerobic bioprocesses scale-up.
In this work, the scale-up procedure is first developed for a general bioreactor, revealing the Hankel matrix singular value decomposition as a key tool in the scale-up problem of bioprocesses.Then, the procedure is applied to a polyhydroxybutyrate (PHB) fed-batch bioreactor allowing the process to be scaled up to a 90 L unit, achieving similar process dynamic behavior and PHB production at both laboratory and industrial scales.Here, a simulation comparison between three empirical scale-up criteria (constant volumetric power consumption, impeller tip speed and Reynolds number) and the proposed procedure is evaluated.

Methods
Consider the following matrix equation that describes the mass distribution within a general bioreactor [11].
where x 2 R n represents the state vector, y 2 R m the output vector consisting of a set of key variables from a design viewpoint, and x in 2 R n the input concentration vector.K 2 R nÂl is the matrix containing stoichiometric coefficients (yields) and r Á ð Þ 2 R l the reaction rate vector.q 2 R n represents the mass transfer from the liquid to the gaseous phase and f 2 R n the mass transfer from the gaseous to the liquid phase.Here, D 2 R þ is the dilution rate, equal to the ratio of the input flow rate (F in ) on the volume fraction (V L ).
From ( 1) and ( 2), process variables and parameters are classified into state variables (x), design variables (z), synthesis parameters (p) and design-variable-dependent functions (w).In this sense, the model can be reformulated as follows: where z 2 R p is the design vector comprising all variables that can be freely varied by the designer, meaning that they can be changed according to the new scale requirements.z is composed of each species input concentration and the bioreactor initial liquid volume.In (3), p 2 R h represents the synthesis parameters that are inherent parameters of the process such as each species yield, inhibition and saturation constants, maintenance coefficients, among others.Once synthesis parameters are established at the current scale, they are considered to be constant for the scale-up.Also, the product Mz 2 R n represents the input concentration vector, where M 2 R nÂp selects the corresponding input concentrations from z. Here, K is purely composed of process synthesis parameters and r depends on the state variables and synthesis parameters.Finally, q; f and D depend on x; z; p and may depend on each other (as shown in the next section), thus they can be considered as designvariable-dependent functions (w).According to this, w depends on z; i.e. each of the entries of w can be written as an explicit function of z: Typically, parameters such as each species flow rate and volumetric mass transfer coefficients belong to this vector.It is worth clarifying that z is vector of fixed parameters that the designer may change according to scale increments, on the contrary, w is a vector of time-varying algebraic expressions that may change throughout scale changes.
Assuming that under usual working conditions this process operates along the operating trajectory x N while is driven by the system input z N and that the nonlinear system operates inside a small neighborhood of the nominal trajectory, the process dynamics, described by ( 3) and ( 4), can be approximated by the first term of its Taylor series expansion.
where A c ; B c and C c are the Jacobian matrices of the system given by ( 7)- (9).
Here, subscript N represents the nominal operating trajectory.Notice that since the bioprocess model is linearized along the operating trajectory, rather than a fixed point, linear A c ; B c and C c matrices are, in general, time varying.B c and C c should be modified, using (10) and (11), to make both design and output variables dimensionless and normalized, guaranteeing for both of them to be homogeneous in magnitude.A c is not altered because the Hankel matrix is a tool that only considers inputs and outputs of the system [3], so any mathematical operation done over x will be annulled during the Hankel matrix calculation [2].
where subscripts max and min are the maximum and minimum values of z j and y i in each case.This type of transformation is called ''scaling'' in [1,44].However, the use of this word is avoided within this work to prevent any confusion when mentioning the words ''scale-up'' associated to scale increments.B c and C c transformation enables avoiding inaccuracy in the model from disparate values and units in process inputs and outputs, which always results in highly large and small entries into both matrices and, hence in an ill-conditioned problem [2,25].
The linear model consisting of ( 5) and ( 6) is conveniently discretized as shown: Here, are the process discrete matrices.In addition, the discrete model composed of ( 12) and ( 13) can be rewritten as shown in ( 14) [2,28,38].
Singular values of the Hankel matrix are closely related to the controllability and observability of the system [51,53].Therefore, SVD of H as shown in (18) provides additional information regarding the controllability and observability of the process, which can be used to quantify the importance of each state in the corresponding input-output system described by H [2,15].
where the matrices U 2 R nmÂnm and V 2 R npÂnp are the column (output) and row (input) spaces of H: Also, the diagonal elements of R 2 R nmÂnp are the singular values (r ii ) of H: From a physical perspective, the matrix of singular values (R) provides the information intensity of the system represented by H; where the highest singular value contains most of the system information [34].This means that if the Hankel singular values decrease rapidly, most of the input-output behavior is provided by the first few states [15,57].Therefore, by means of singular values analysis, the information contained in H can be reconstructed just by representing the system with the input-output variables related to the highest singular values [2,53].
Given that SVD of the Hankel matrix determines in a qualitative way whether a process is controllable/observable or not, as a way to quantify these process information extracted from (18), the output impactability index of each output variable (OII y k ) was introduced in [2].Considering that U represents the output space of H as can be noticed from (18), and hence each column of U is related to one output (left singular vector) [34,53], OII y k is defined as (19).The rationale behind this index is based on the concept of Euclidean norm, therefore, the sum of each squared output singular vector entries related to each output at each time instant U 2 kþms;i ðjÞ weighted by the corresponding squared singular value r 2 ii ðjÞ determines the importance of each output variable in the process [2].
ii ðjÞ where k ¼ 1; 2; . ..; m and r is the rank of H, i.e. the number of non-zero r ii .It is the system order and defines the dimension of the controllable and observable subspace [56].OII y k represents the impactability of the process design variables (z) as a whole over a kth given output variable (y k ).Here, the most impacted output variable (main dynamics) is the y k with the highest OII and it corresponds to the process governing dynamics.Therefore, based on the OII calculation, a quantitative hierarchical relation among the output variables of the process can be established.This relation is here called dynamics hierarchy and it determines the degree of importance of each dynamics along the process.Besides, based on the determination of the dynamics hierarchy, the point of the fed-batch trajectory where the process should be scaled up is established, called here the critical point of the operating trajectory.This point is determined from the OII curve for the main dynamics (governing dynamics) as the maximum value of the index for such dynamics.Given that the main dynamics is the output variable with the highest OII, slight changes in such dynamics may affect considerably the process evolution and, therefore, the point where the OII reaches its maximum value represents the time when the process output variables are most affected by the design variables, meaning that this is the best point of the operating trajectory for scaling up the process.
On the other hand, the critical point also represents the time instant of the operating trajectory where specific process requirements are maximum, e.g.where the mass transfer is governing the overall process progress.This is because given that z is the vector comprising all variables that can be freely varied by the designer during the scale-up and the main dynamics corresponds to the governing process dynamics, the critical point corresponds to the time when the process requirements are maximum according to the design variables current values.So, for any point with an OII less than the OII computed at the critical point, the process requirements are fulfilled, considering that those requirements were fulfilled at the critical point.
The idea underlying this proposal is to establish a new unit design in agreement with the process dynamics.To do so, the OII at the critical point must remain constant as the scale is increased and each design-variable-dependent function (w) should be given by a valid equation at both scales, i.e. an expression that maintains the OII of the critical point invariant throughout the scale-up.Bearing this in mind, based on this methodology implementation, it can be evaluated if certain scale-up criterion will deliver a new scale design with similar dynamics to the current scale.Therefore, if the OII value at the critical point remains constant throughout the scale-up, the current scale dynamic behavior is transferred to the new scale.As a result, the potential of this index can be exploited to improve the scale-up problem of biotechnological processes as shown in the next section.

Results and discussion
The process of polyhydroxybutyrate (PHB) production comprises two stages [6,19,32]: (i) biomass growth and (ii) PHB production, as can be seen in Fig. 1.
During the initial stage, required nutrients that enable the biomass growth (carbon, nitrogen and oxygen sources) are supplied, allowing biomass concentration to rise up to the desired level [16].In the second stage, the nitrogen source input is stopped, disrupting biomass growth, allowing the excess carbon to be used in the PHB production [32].Here, it is assumed that the kinetic mechanism, described by ( 20)-( 22), governs the PHB synthesis [33].
where (20) describes the biomass growth on glucose, (21) the biomass growth on PHB and ( 22) the PHB production.In ( 20)-( 22), substrate, nitrogen, oxygen, biomass, PHB and carbon dioxide are represented by S, N, O, X, P and CO 2 , respectively.Yield coefficients involved in ( 20)-( 22) together with all variables and parameters within the model are defined in the Nomenclature section.The PHB production model, ( 23)-( 45), is borrowed from [33], but including the oxygen dynamics at both liquid and gaseous phases [8,40].From ( 23)- (29), it can also be noticed that the mathematical structure of the proposed model belongs to the general matrix model of ( 1) and (2).
Fig. 1 Fed-batch bioreactor process Substrate, nitrogen and air flow rates are calculated using ( 43)-( 45) respectively, where ( 43) and ( 44) are open loop control laws that maintain both substrate and nitrogen concentrations constant.A practical application of this control strategy can be seen in [32] On the other hand, ( 45) is a closed loop law that regulates the dissolved oxygen concentration at O L;sp , where eðtÞ ¼ O L;sp À O L considering that O L;sp ¼ 0:55O H L for the first stage and O L;sp ¼ 0:30O H L for the second stage.Figure 2 shows the dynamic evolution of the most important process variables compared with the experimental data reported in [33], where X t ¼ X þ P and the subscript exp represents the experimental data.
Here, it can be noticed that the proposed model represents the experimental data, as was reported in [33].Values of model parameters are reported in [33], the additional ones related to the oxygen dynamics are listed in Table 1.
Notice that since O G;F is the oxygen concentration in the air feeding stream to the bioreactor, such concentration could be increased either mixing the air stream with pure oxygen or feeding pure oxygen to the bioreactor.
The process model is linearized along the operating trajectory as shown in ( 5) and (6), where A c , B c and C c are given by ( 7)- (8).Output and input variables are normalized using (10) and (11), considering that y i;min and y i;max are minimum and maximum values of each y i along the process and, z j;min and z j;max are AE10 % of their nominal values, since minor changes are expected for these limits while the scale is increased because the process dynamic behavior is transferred from the current to the new scale using the proposed methodology [30].Then, the model is discretized as shown in ( 12) and ( 13) and, O, C and H are computed using ( 15)- (17).Subsequently, H is factorized via SVD using (18) to, finally, compute the OII y k using (19), where rankðHÞ ¼ r ¼ 5. Here, the output variables are considered to be the process state variables in order to determine which of them is the most important from a design viewpoint.It is worth highlighting that the SVD of H allows the designer to determine the effect of the design K O 0:000118 g L [42] m O 0:008987 N p 5 dimensionless [12] O G;F 0:269 g L [45] Fig. 2 State variables profiles at the current scale (input) variables over each state (output) variable.In this sense, Fig. 3 shows the OII y k profiles, where it can be seen that O L is the most impacted dynamics by the design variables of the process (variable with highest index) and, that the critical point of the batch is located at the end of the first process stage, where the oxygen requirement is maximum for the process.This point corresponds to OII O L ¼ 0:019 at t ¼ 31 h.Therefore, the process should be scaled up at this point using a valid expression for each design-variable-dependent function of w.Here, F S , F N and F O are computed with ( 43)- (45).On the other hand, K L a is computed using four criteria for the sake of comparison: (1) process requirements (herein proposed), (2) volumetric power consumption, (3) impeller tip speed and, (4) Reynolds number.
According to the first criterion, K L a should be established by means of the required amount of oxygen for biomass growth.Therefore, (50) fulfills this condition because this expression is based on the mass balance for the oxygen in the gaseous phase, considering the oxygen demand at the liquid phase.
Notice that K L a depends on F O .The second criterion determines the new scale unit design holding the volumetric power consumption from current to new scale [17,47].In this sense, ( 51) is used to calculate the power requirement and (39) to establish the volumetric mass transfer coefficient at the new scale [26,54].
In the third criterion, it is stated that the impeller tip speed (V t ) must remain constant as the scale is increased.
Considering that V t ¼ pN i D i [41,55], ( 52) is used to determine the impeller speed and (39) to compute K L a at the new scale [17].
The last criterion establishes that the process should be scaled up maintaining the same Reynolds number (Re) [17,26].Taking into account that Re is given by ( 53), this empiric rule can be simplified into (54).Here, K L a is also computed with (39).
Equations ( 50)-( 54) are evaluated at the critical point.It must be noticed that (50) does not fix the reactor dimensions in opposition to (39) that provides a specific geometry of the process unit.Therefore, reactor dimensions are similar to the current scale design when using empirical rules, meaning that K L a is fixed by the unit geometry while when using (50) the process unit can be completely redesigned at the new scale in agreement with the oxygen requirements (process requirements).Bearing this in mind, Table 2 summarizes the volumetric oxygen transfer coefficient, reactor dimensions, and the OII O L for all the previously described cases at the process critical point.Here, it is worth highlighting that K L a and OII O L are equal when the process is scaled up maintaining the process requirements, supporting the widely used criterion of keeping K L a constant [7,10,17,22,47,52].It can also be noticed from Table 2 a larger K L a value was estimated when the volumetric power consumption is kept constant.This means that a larger process unit than the required one was designed using this criterion.For the other two cases (constant V t and Re), on the other hand, a smaller K L a value was established at the new scale, meaning that a smaller process unit than the required one was designed.This deduction can also be established from Figs. 4 and 5, where a comparison of the dissolved oxygen dynamics and the air flow rate is done for all cases.
From Fig. 4, it can be seen that curves for the process requirements and constant power overlap the current scale curve.For the former, it is demonstrated that the process requirements are identical at both scales (see K L a in Fig. 3 Output impactability index at the current scale Table 2), revealing that the oxygen dynamics governs the overall rate of the process at both scales.However, although the oxygen dynamics is similar for the latter case (constant power), Fig. 5 shows that the required air flow rate is less than the established at the current scale, confirming that the unit was oversized (see K L a in Table 2).
It can also be noticed in Fig. 5 that when scaling up the process maintaining the impeller tip speed, the air flow rate reaches its maximum value during the last period of the first process stage, deteriorating the dissolved oxygen behavior (see Fig. 4) and corroborating that this unit was undersized (see K L a in Table 2).During this time, any disturbance introduced to the process cannot be countered by the controller.Here, the same controller parameters were used for all cases.
On the other hand, regarding the case of constant Reynolds number, the air flow rate reaches its maximum value from the beginning of the process, meaning also that this unit was undersized (see K L a in Table 2).This process behavior demonstrates that maintaining the same Reynolds number as scale-up criterion hardly ever results in an adequate unit design as was stated in [26], since by means of this rule the estimated gassed input power is always lower than the effectively required.
In addition, a comparison of the ratio of oxygen concentration in the gas phase (O g ) on PHB concentration (P) for all cases is done in Fig. 6.
It can be noticed that the new scale profile for the constant power case is highly different from the current scale curve.This behavior ratifies that the unit was oversized.Here, given that this ratio is smaller at the new scale than the current one, it is possible to conclude that K L a was overestimated when using this criterion.In opposition, O g P is greater for the other two cases (constant tip speed and  Reynolds number) than the current scale profile, confirming that both units were undersized and that K L a was underestimated when using these criteria.Finally, from Fig. 6, it can be seen that the profile for the process requirements case overlaps the current scale curve, revealing that the process has the same dynamic behavior at new scale when the process is scaled up using the proposed scale-up approach (see also Figs. 4, 5).Here, the process reaches the same ratio of oxygen concentration in the gas phase on PHB concentration from the current scale at new scale, which demonstrates that maintaining the same value of the volumetric oxygen transfer coefficient when increasing the process scale is the best method for scaling up this type of process.

Conclusions
This work presents a general approach to the scale-up of aerobic fed-batch bioprocesses, based on model of the process dynamics.The methodology herein proposed is based on the singular value decomposition of the process Hankel matrix taking advantage of the relation of Hankel singular values with the process controllability and observability.By means of this proposal, the designer is able to determine the most relevant process variables through the calculation of the Output Impactability Index.This index allows the designer to determine the point where the process should be scaled up to fulfill its dynamic requirements.Also, by means of the Output Impactability Index calculation, it is possible to establish if two or more designed units can carry out the same process with the same performance targets.
A fed-batch fermenter was scaled up from 3 to 90 L in a realistic simulation environment.As a result, the scale factors for keeping the same ratio of oxygen concentration in the gas phase (O g ) on polyhydroxybutyrate (PHB) concentration, with the same dissolved oxygen dynamics, at both scales were found.Here, the most impacted dynamics by the design variables is the dissolved oxygen concentration (O L ) confirming that the oxygen dynamics governs the overall rate of the bioprocess.In addition, it was established that the volumetric oxygen transfer coefficient remains constant through the scale-up, which also ratifies that the oxygen transport from the gas/liquid interface to the bulk of the liquid is a key phenomenon when scaling-up this kind of fermentation.

Fig. 4 Fig. 5
Fig. 4 Comparison of the dissolved oxygen concentrations at the current and new scale

Table 1
Additional model parameters

Table 2
Comparison of the new and current scale unit design for all criteria