Model-Based Evolutionary Operation Design for Batch and Fed- Batch Antibiotic Production Bioprocesses

Abstract

The control policy determination for batch and fed-batch antibiotic production bioprocesses is an important practical issue due to the high added value of these bioproducts. Since it is highly desirable to optimize the antibiotic production, several methods have been proposed aimed at this objective. Once having a mathematical model for the bioprocess, the optimization problem can be formulated within the framework of Pontryagin's maximum principle and of the optimal control theory to determinate the best control trajectory for certain key manipulated variables, such as temperature, pH, and substrate feed rate. In this chapter, applications of these model-based techniques to optimize and control antibiotics production bioprocesses are reviewed and new aspects are emphasized. The cases analyzed included the optimization of the substrate feed rate in a fed-batch reactor and of the temperature in a batch reactor during penicillin fermentations. The main contributions of this study were: (i) the proposition of a different procedure for calculating the second switching time of substrate feed rate, (ii) the application of simpler numerical methods to solve the two-point boundary-value problem associated with the temperature profile optimization, and (iii) the demonstration that the non-isothermal operation is more productive in antibiotic than the operation under constant temperature.

1. Introduction

Improvement in the productivity of many submerged fermentation processes is carried out by manipulating nutritional and physical parameters such as medium composition, agitation speed, aeration rate, pH, and temperature [1, 2]. Although the attainment of optimal conditions for a multivariable fermentation process is often tedious, it is possible to undertake a rational procedure by using statistical experimental designs [2].

Experimental designs can be divided in two distinct groups [3, 4, 5, 6]: (i) model-based experimental designs and (ii) statistical experimental designs. In model-based experimental designs, predictions of a mathematical model are used to determine how an experiment or process should be performed, whereas, with statistical experimental designs, these model predictions are not explicitly required.

The optimization and operation of fermentation processes play a key role in the biotechnology industry due to heavy competition among companies. Secondary metabolites, such as antibiotics and other pharmaceutical products, represent an important added value; therefore, improvements in the production of these bioproducts are of great interest to industries. To achieve high-performance operations, the optimization of manipulated variables that affect the fermentation process becomes a significant task.

In general, optimization problems can be classified in two categories: set-point and profile optimizations [7]. Set-point optimization problems involve finding the best set of values of manipulated variables that lead to the maximization of performance indexes [7]. Profile optimization consists of determining temporal or spatial functions (profiles), rather than a point in n-dimensional space, which lend an optimal value to the performance index [7].

In antibiotic fermentation, it is well known that the temperature and pH for the maximum rate of antibiotic production are different from those for the maximum rate of cell growth [8]. In this sense, the implementation of temperature and pH profiles plays an important role in significant improvements in antibiotic production bioprocesses [8].

Since the primary goal of a fermentation process is the cost-effective production of bioproducts, it is important to select the more appropriate operating mode that allows the production of the desired product at a high concentration with a high productivity and yield [9]. Fed-batch bioprocesses have been widely employed for the production of various bioproducts, including primary and secondary metabolites [9]. In the particular case of secondary metabolites, such as antibiotics, the interaction between growth metabolism and product biosynthesis is critically affected by growth-limiting nutrient concentrations. Since both the underfeeding and the overfeeding of nutrients are detrimental to cell growth and product formation, due to the occurrence of phenomena such as cell starvation and catabolite repression, establishing a suitable feeding strategy is crucial in fed-batch bioprocesses [9, 10].

A particular time sequence of control variables may be required in order to conduct the bioprocess over time in a trajectory that provides the greatest productivity. This can lead to complex optimal time profiles for the control variables, which are sometimes impossible to be determined purely experimentally. Thus, appropriate mathematical and numerical methods can be applied for the determination of these profiles in order to reduce the experimental effort and the required time for optimization.

The search for the optimal pH, temperature, and substrate feed-rate profiles in batch and fed-batch antibiotic fermentation is a typical problem of optimization and evolutionary operations for which the use of kinetic models and powerful mathematical techniques is essential for their solution [7, 10, 11]. According to Rani and Rao [12], several approaches for the determination of optimal time profiles for control variables have been reported in the literature [13, 14, 15]. In these reports, the optimization problem is generally formulated on the basis of Pontryagin’s maximum principle, taking as a starting point a phenomenological mathematical model of the bioprocess. For simple mathematical models, the problem can be solved analytically, from the Hamiltonian of the system, by applying an iterative scheme on the control variable to determine the optimal control profile [8, 16, 17, 18, 19].

In this chapter, two studies on the optimization and the evolutionary operation of antibiotic production bioprocesses are revisited, and new results are obtained and highlighted. Such studies report mathematical models of bioprocesses, in conjunction with Pontryagin’s maximum principle, to optimize the substrate feed-rate profile for a fed-batch bioreactor and the temperature profile for a batch fermentation in order to maximize the production of antibiotic. The fundamentals of Pontryagin’s maximum principle, when applied to the cases analyzed, are also presented.

The aim is to provide a theoretical basis for the application of a model-based methodology that can be used for the optimization and control of bioprocesses from other antibiotics and secondary metabolites with a broad structural diversity and therapeutic activity, including antibacterial, antifungal, antiviral, antitumor, immunosuppressive, antihypertensive, and antihypercholesterolemic compounds.

2. Case studies: batch and fed-batch antibiotic production bioprocesses

During batch and fed-batch bioprocesses, the state variables (cell, substrate, oxygen and product concentrations, temperature, and pH) change significantly, from initial to final values. This dynamic behavior motivates the development of optimization methods to find the optimal time trajectories for the control variables in order to improve the performance of these bioprocesses.

Two case studies on the optimization of control variables in batch and fed-batch antibiotic production bioprocesses are revisited, and additional results are obtained and presented. The cases studied are those reported by Costa [20] and Constantinides and Mostouffi [17], concerning the optimization of the substrate feed rate in a fed-batch reactor and the temperature in a batch reactor, respectively. These cases are presented and detailed in the following sections.

2.1. Case study #1: determination of the optimal substrate feed-rate profile in fed-batch bioreactor for penicillin production

In this case study, the bioprocess of penicillin production by Penicillium chrysogenum is described by the mathematical model presented by Costa [20], which is based on the classical model proposed by Bajpai-Reuss for penicillin fermentation. In this model, the specific growth rate (μ) takes into account diffusional limitations that occur in the filamentous fungal biomass, as described by the Contois model. The specific rate of product formation (π) considers that penicillin production is repressed by high substrate concentrations (catabolic repression), being modeled by the Andrews equation. Penicillin degradation by hydrolysis is also considered, assuming first-order kinetics for this reaction. The specific rate of substrate consumption (σ) is represented by the Herbert-Pirt generalized model, whereby the substrate is consumed for cell growth and maintenance and for product formation. The equations of the full mathematical model are as follows (note that if the dilution rate is used as the control variable, the total mass-balance equation is not required for the optimization problem formulation, thus reducing the equation system dimension):

$\frac{dX}{dt}=\mu X-Xu;\mu =\frac{{\mu }_{m}S}{BX+S}$ E1-2

$\frac{dS}{dt}=-\sigma X+(S_f-S)u;\sigma =\frac{\mu }{{\text{Y}}_{\text{X/S}}}+\frac{\pi }{Y_{P/S}}+m$ E3-4

$\frac{dP}{dt}=(\pi X-k_{h}P)-Pu;\pi =\frac{{\pi }_{m}S}{k_m+S+\frac{S^2}{k_i}}$ E5-6

where

• X, S, and P denote the concentrations of cell, substrate, and product, respectively;

• μ, σ, and π are the specific rates of cell growth, substrate consumption, and product formation;

• μ_m, B, Y_X/S, Y_P/S, m, π_m, k_m, and k_i are the parameters of the mathematical model, including kinetic and yield parameters;

• $u=D=F/V:u=\text{control variable};D=\text{dilution rate};F=\text{feed rate};V=\text{culture volume}$ E7

In matrix notation:

$\frac{d{\displaystyle \underset{¯}{X}}}{dt}={\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)u$ E8

where

${\displaystyle \underset{¯}{X}}=\left[\begin{array}{c}X\\ S\\ P\end{array}\right];{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)=\left[\begin{array}{c}\mu X\\ -\sigma X\\ \pi X-k_{h}P\end{array}\right];{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)=\left[\begin{array}{c}-X\\ (S_f-S)\\ -P\end{array}\right]$ E9-11

The constraints imposed on the control variable u are as follows:

$u_{\min }\le u\le u_{\max }$ E12

where u_min = 0 and u_max = F_max / V_.

Another constraint concerns the maximum volume of culture (final volume), i.e., V(t_f) = V_max = V_f, where t_f is the final processing time.

The initial conditions are given by

$X(0)=X_0\text{;}S(0)=S_0;P(0)=P_0;V(0)=V_0$ E13

The objective of the optimization/control problem is to determine the optimal time profile for the control variable that maximizes the antibiotic concentration at the end of the bioprocess.

According to the Pontryagin’s maximum principle, the optimal profile must maximize the Hamiltonian, given by

$H={{\displaystyle \underset{¯}{\lambda }}}^T[{\displaystyle \underset{¯}{f}}({\displaystyle \underset{¯}{X}})+{\displaystyle \underset{¯}{g}}({\displaystyle \underset{¯}{X}})u];{{\displaystyle \underset{¯}{\lambda }}}^T=[{\lambda }_1{\lambda }_2{\lambda }_3]$ E14

or

$H=H_0\left(t\right)+\phi \left(t\right)u\{\begin{array}{l}{H}_0\left(t\right)={{\displaystyle \underset{¯}{\lambda }}}^T{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)\\ \phi \left(t\right)={{\displaystyle \underset{¯}{\lambda }}}^T{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\end{array}$ E15

Since ${\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)=\left[\begin{array}{l}{f}_1\\ f_2\\ f_3\end{array}\right]$ and ${\displaystyle \underset{¯}{\text{g}}}\left({\displaystyle \underset{¯}{X}}\right)=\left[\begin{array}{l}{g}_1\\ g_2\\ g_3\end{array}\right]$ , the following equations are obtained:

$H_0\left(t\right)={\lambda }_1{f}_1+{\lambda }_2{f}_2+{\lambda }_3{f}_3={\lambda }_1\mu X-{\lambda }_2\sigma X+{\lambda }_3\left(\pi X-kP\right)$ E16

$\varphi \left(t\right)={\lambda }_1{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3=-{\lambda }_{1}X+{\lambda }_2\left(S_f-S\right)-{\lambda }_{3}P$ E17

For optimal control, it is established that

• If $\phi \left(t\right)>0$ , then $u=u_{\max }$ .

• If $\phi \left(t\right)<0$ , then $u=u_{\min }$ .

• If $\phi \left(t\right)=0,$ then $u=u_{\sin }$ .

$\text{In the singular interval: }H(t)=0;\varphi (t)=0;H_0(t)=0$ E18-20

Since ${\displaystyle \underset{¯}{\lambda }}={\displaystyle \underset{¯}{\lambda }}\left(t\right)$ , the following equations can be developed:

$\frac{d}{dt}{\displaystyle \underset{¯}{\lambda }}=-\frac{\partial H}{\partial {\displaystyle \underset{¯}{X}}}=-\frac{\partial }{\partial {\displaystyle \underset{¯}{X}}}\text{[}{{\displaystyle \underset{¯}{\lambda }}}^T(f(x)+g(x)u)\text{]}=-{{\displaystyle \underset{¯}{\lambda }}}^T{\displaystyle \underset{{\left(\frac{d}{dt}{\displaystyle \underset{¯}{\lambda }}\right)}^T}{{\displaystyle \underbrace{\left(\frac{d}{d{\displaystyle \underset{¯}{X}}}{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+u\frac{d}{d{\displaystyle \underset{¯}{X}}}{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right)}}}}$ E21

$\varphi \left(t\right)={{\displaystyle \underset{¯}{\lambda }}}^{T}g\left({\displaystyle \underset{¯}{X}}\right)\Rightarrow \frac{d\varphi }{dt}={\left(\frac{d}{dt}{\displaystyle \underset{¯}{\lambda }}\right)}^T{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)+{{\displaystyle \underset{¯}{\lambda }}}^T\frac{d}{dt}{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)$ E22

$\frac{d}{dt}{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)=\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right)\frac{d{\displaystyle \underset{¯}{X}}}{dt}=\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right)\text{[}{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)u\text{]}$ E23

$\frac{d}{dt}{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)=\frac{d}{d\left({\displaystyle \underset{¯}{X}}\right)}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)u$ E24

Substituting Eqs. (21) and (24) into Eq. (22) gives

$\frac{d\phi }{dt}=-{{\displaystyle \underset{¯}{\lambda }}}^T\left(\frac{d}{d{\displaystyle \underset{¯}{X}}}{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+u\frac{d}{d{\displaystyle \underset{¯}{X}}}{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)+{{\displaystyle \underset{¯}{\lambda }}}^T\left(\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)u\right)$ E25

$\frac{d\phi }{dt}={{\displaystyle \underset{¯}{\lambda }}}^T\left(\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)-\frac{d}{d{\displaystyle \underset{¯}{X}}}\left({\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)\right){\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)\right).$ E26

By developing the matrices indicated in the previous equation, one obtains

$\frac{d\phi }{dt}=\text{[}{\lambda }_1{\lambda }_2{\lambda }_3\text{]}*\text{[}\begin{array}{l}\frac{\partial g_1}{\partial X_1}\frac{\partial g_1}{\partial X_2}\frac{\partial g_1}{\partial X_3}\\ \frac{\partial g_2}{\partial X_1}\frac{\partial g_2}{\partial X_2}\frac{\partial g_2}{\partial X_3}\\ \frac{\partial g_3}{\partial X_1}\frac{\partial g_3}{\partial X_2}\frac{\partial g_3}{\partial X_3}\end{array}\text{]}\text{ *}\left[\begin{array}{l}{f}_1\\ f_2\\ f_3\end{array}\right]\text{- }\text{[}\begin{array}{l}\frac{\partial f_1}{\partial X_1}\frac{\partial f_1}{\partial X_2}\frac{\partial f_1}{\partial X_3}\\ \frac{\partial f_2}{\partial X_1}\frac{\partial f_2}{\partial X_2}\frac{\partial f_2}{\partial X_3}\\ \frac{\partial f_3}{\partial X_1}\frac{\partial f_3}{\partial X_2}\frac{\partial f_3}{\partial X_3}\end{array}\text{]}\text{ *}\left[\begin{array}{l}{g}_1\\ g_2\\ g_3\end{array}\right]$ E27

where

$g_1=-X;g_2=(S_f-S);g_3=-P;X_1=X;X_2=S;X_3=P$ E28-33

$\frac{\partial g_1}{\partial X_1}=-1;\frac{\partial g_1}{\partial X_2}=0;\frac{\partial g_1}{\partial X_3}=0$ E34-36

$\frac{\partial g_2}{\partial X_1}=0;\frac{\partial g_2}{\partial X_2}=-1;\frac{\partial g_2}{\partial X_3}=0$ E37-39

$\frac{\partial g_3}{\partial X_1}=0;\frac{\partial g_3}{\partial X_2}=0;\frac{\partial g_3}{\partial X_3}=-1$ E40-42

$f_1=\mu X;f_2=-\sigma X;f_3=\pi X-k_{h}P;\mu =\mu \left(S,X\right);\sigma =\sigma \left(S,X\right);\pi =\pi \left(S\right)$ E43-48

$\frac{\partial f_1}{\partial X_1}=\mu +X{{\mu }^{\prime }}_X;\frac{\partial f_1}{\partial X_2}=X{{\mu }^{\prime }}_S;\frac{\partial f_1}{\partial X_3}=0$ E49-51

$\frac{\partial f_2}{\partial X_1}=-\left(\sigma +X\sigma {\text{'}}_X\right);\frac{\partial f_2}{\partial X_2}=-X{{\sigma }^{\prime }}_S;\frac{\partial f_2}{\partial X_3}=0$ E52-54

$\frac{\partial f_3}{\partial X_1}=\pi ;\frac{\partial f_3}{\partial X_2}=X{{\pi }^{\prime }}_S;\frac{\partial f_3}{\partial X_3}=-k_h$ E55-57

Then

$\frac{d\varphi }{dt}=\text{[}{\lambda }_1{\lambda }_2{\lambda }_3\text{]}*\left[\begin{array}{ccc}-1& 0& 0\\ 0& -1& 0\\ 0& 0& -1\end{array}\right]*\left[\begin{array}{c}{f}_1\\ f_2\\ f_3\end{array}\right]-\left[\begin{array}{ccc}\mu +X{{\mu }^{\prime }}_X& X{{\mu }^{\prime }}_S& 0\\ -\left(\sigma +X{{\sigma }^{\prime }}_X\right)& -X{{\sigma }^{\prime }}_S& 0\\ \pi & X{{\pi }^{\prime }}_S& -k_h\end{array}\right]*\left[\begin{array}{c}{g}_1\\ g_2\\ g_3\end{array}\right]$ E58

$\frac{d\varphi }{dt}=\text{[}{\lambda }_1{\lambda }_2{\lambda }_3\text{]}*\left[\begin{array}{c}-f_1\\ -f_2\\ -f_3\end{array}\right]-\left[\begin{array}{c}{g}_1\left(\mu +X{{\mu }^{\prime }}_X\right)+g_{2}X{{\mu }^{\prime }}_S\\ -g_1\left(\sigma +X{{\sigma }^{\prime }}_X\right)-g_{2}X{{\sigma }^{\prime }}_S\\ g_1\pi +g_{2}X{{\pi }^{\prime }}_S-g_3{k}_h\end{array}\right]$ E59

$\frac{d\varphi }{dt}={\lambda }_1\left(-f_1-g_1\left(\mu +X{{\mu }^{\prime }}_X\right)-g_{2}X{{\mu }^{\prime }}_S\right)+{\lambda }_2\left(-f_2+g_1\left(\sigma +X{{\sigma }^{\prime }}_X\right)+g_{2}X{{\sigma }^{\prime }}_S\right)+{\lambda }_3\left(-f_3-g_1\pi -g_{2}X{{\pi }^{\prime }}_S+g_3{k}_h\right)$ E60

$\frac{d\varphi }{dt}={\lambda }_1\left(-\mu X+X\left(\mu +X{{\mu }^{\prime }}_X\right)+\left(S_f-S\right)X{{\mu }^{\prime }}_S\right)+{\lambda }_2\left(\sigma X-X\left(\sigma +X{{\sigma }^{\prime }}_X\right)+\left(S_f-S\right)X{{\sigma }^{\prime }}_S\right)+{\lambda }_3\left(-\left(\pi X-k_{h}P\right)+\pi X-\left(S_f-S\right)X{{\pi }^{\prime }}_S-k_{h}P\right)$ E61

$\frac{d\varphi }{dt}={\lambda }_1\left({{\mu }^{\prime }}_X{X}^2-\left(S_f-S\right)X{{\mu }^{\prime }}_S\right)+{\lambda }_2\left(-{{\sigma }^{\prime }}_X{X}^2+\left(S_f-S\right)X{{\sigma }^{\prime }}_S\right)+{\lambda }_3\left(-\left(S_f-S\right)X{{\pi }^{\prime }}_S\right)$ E62

In the singular interval:

$\frac{d\varphi }{dt}=0$ E63

$H_0\left(t\right)={\lambda }_1{f}_1+{\lambda }_2{f}_2+{\lambda }_3{f}_3=0$ E64

$\varphi \left(t\right)={\lambda }_1{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3=0$ E65

From Eq. (64):

${\lambda }_1{f}_1+{\lambda }_2{f}_2+{\lambda }_3{f}_3=0\Rightarrow {\lambda }_1=-\frac{\left({\lambda }_2{f}_2+{\lambda }_3{f}_3\right)}{f_1}$ E66

Substituting Eq. (66) into Eq. (65) results in the following equation:

$\begin{array}{l}-\frac{\left({\lambda }_2{f}_2+{\lambda }_{3}f\right)}{f_1}{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3=0\Rightarrow -\frac{{\lambda }_2{f}_2}{f_1}{g}_1-\frac{{\lambda }_3{f}_3}{f_1}{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3=0\Rightarrow \\ \left(g_2-\frac{f_2}{f_1}{g}_1\right){\lambda }_2=\left(\frac{f_3}{f_1}{g}_1-g_3\right){\lambda }_3\Rightarrow {\lambda }_2=\frac{\left(\frac{f_3}{f_1}{g}_1-g_3\right){\lambda }_3}{\left(g_2-\frac{f_2}{f_1}{g}_1\right)}\Rightarrow {\lambda }_2=\underset{\beta }{\underbrace{\left(\frac{f_3{g}_1-f_1{g}_3}{f_1{g}_2-f_2{g}_1}\right)}}{\lambda }_3\end{array}$ E67

By introducing the expression of λ₂ into the expression of λ₁, one obtains

${\lambda }_1={\displaystyle \underset{\alpha }{{\displaystyle \underbrace{-\text{[}\frac{\left(\frac{f_3{g}_1-f_1{g}_3}{f_1{g}_2-f_2{g}_1}\right)f_2+f_3}{f_1}\text{]}}}}}{\lambda }_3$ E68

Substituting λ₁=αλ₃ and λ₂=βλ₃ into the expression of $d\varphi dt$ in the singular interval provides

$\frac{d\varphi }{dt}={\lambda }_3\alpha \left({{\mu }^{\prime }}_X{X}^2-\left(S_f-S\right)X{{\mu }^{\prime }}_S\right)+{\lambda }_3\beta \left(-{{\sigma }^{\prime }}_X{X}^2+\left(S_f-S\right)X{{\sigma }^{\prime }}_S\right)+{\lambda }_3\left(-(S_f-S)X{{\pi }^{\prime }}_S\right)=0$ E69

$\frac{d\varphi }{dt}={\lambda }_3{\displaystyle \underset{Q\left({\displaystyle \underset{¯}{X}}\right)}{{\displaystyle \underbrace{\text{[}\alpha \left({{\mu }^{\prime }}_X{X}^2-\left(S_f-S\right)X{{\mu }^{\prime }}_S\right)+\beta \left(-{{\sigma }^{\prime }}_X{X}^2+\left(S_f-S\right)X{{\sigma }^{\prime }}_S\right)+\left(-\left(S_f-S\right)X{{\pi }^{\prime }}_S\right)\text{]}}}}}=0$ E70

$\frac{d\varphi }{dt}={\lambda }_{3}Q\left({\displaystyle \underset{¯}{X}}\right)=0\Rightarrow Q\left({\displaystyle \underset{¯}{X}}\right)=0$ E71

The equation $Q\left({\displaystyle \underset{¯}{X}}\right)=0$ is the expression of the singular arc, which is independent of the adjoint variables λ₁, λ₂, λ₃. In the expression of $Q\left({\displaystyle \underset{¯}{X}}\right)$ , the indicated derivatives are given by

$\mu =\frac{{\mu }_{m}S}{BX+S}\Rightarrow {{\mu }^{\prime }}_X=-\frac{{\mu }_{m}BS}{{\left(BX+S\right)}^2};{{\mu }^{\prime }}_S=\frac{\left(BX+S\right){\mu }_m-{\mu }_{m}S}{{\left(BX+S\right)}^2}$ E72-73

$\pi =\frac{{\pi }_{m}S}{k_m+S+\frac{S^2}{k_i}}\Rightarrow {{\pi }^{\prime }}_S=\frac{\left(k_m+S+\frac{S^2}{k_i}\right){\pi }_m-{\pi }_{m}S\left(1+\frac{2}{k_i}S\right)}{{\left(k_m+S+\frac{S^2}{k_i}\right)}^2}$ E74

$\sigma =\frac{\mu }{Y_{X/S}}+\frac{\pi }{Y_{P/S}}+m\Rightarrow {{\sigma }^{\prime }}_X=\frac{{{\mu }^{\prime }}_X}{Y_{X/S}};{{\sigma }^{\prime }}_S=\frac{{{\mu }^{\prime }}_S}{Y_{X/S}}+\frac{{{\pi }^{\prime }}_S}{Y_{P/S}}$ E75-76

For the determination of the singular dilution rate (u_sin), one starts from the first derivative of ϕ as follows:

$\frac{d\varphi }{dt}={\lambda }_{3}Q\left({\displaystyle \underset{¯}{X}}\right)\Rightarrow \frac{d^2\varphi }{d{t}^2}={\lambda }_3\frac{d}{dt}Q\left({\displaystyle \underset{¯}{X}}\right)+Q\left({\displaystyle \underset{¯}{X}}\right)\frac{d{\lambda }_3}{dt}$ E77

$\frac{d}{dt}{\displaystyle \underset{¯}{\lambda }}=-\frac{\partial H}{\partial {\displaystyle \underset{¯}{X}}}\Rightarrow \left[\begin{array}{c}d{\lambda }_1/dt\\ d{\lambda }_2/dt\\ d{\lambda }_3/dt\end{array}\right]=-\left[\begin{array}{c}\partial H/\partial X_1\\ \partial H/\partial X_2\\ \partial H/\partial X_3\end{array}\right]$ E78

where

$H=H_0\left(t\right)+\varphi \left(t\right)u={\lambda }_1{f}_1+{\lambda }_2{f}_2+{\lambda }_3{f}_3+\left({\lambda }_1{g}_1+{\lambda }_2{g}_2+{\lambda }_3{g}_3\right)u$ E79

Thus:

$\frac{d{\lambda }_3}{dt}=-\frac{\partial H}{\partial X_3}=-{\lambda }_3\frac{\partial f_3}{\partial X_3}-{\lambda }_3\frac{\partial g_3}{\partial X_3}u$ E80

Substituting the expression of dλ₃/dt into the expression of the second derivative of φ gives

$\frac{d^2\varphi }{d{t}^2}={\lambda }_3\frac{d}{dt}Q\left({\displaystyle \underset{¯}{X}}\right)+Q\left({\displaystyle \underset{¯}{X}}\right)\left(-{\lambda }_3\frac{\partial f_3}{\partial X_3}-{\lambda }_3\frac{\partial g_3}{\partial X_3}u\right)$ E81

When u = u_sin, the second derivative of ϕ is zero, i.e., $d^2\varphi d{t}^2=0$ . Thus,

${\lambda }_3\frac{d}{dt}Q\left({\displaystyle \underset{¯}{X}}\right)+Q\left({\displaystyle \underset{¯}{X}}\right)\left(-{\lambda }_3\frac{\partial f_3}{\partial X_3}-{\lambda }_3\frac{\partial g_3}{\partial X_3}u\right)=0\Rightarrow \frac{d}{dt}Q\left({\displaystyle \underset{¯}{X}}\right)-Q\left({\displaystyle \underset{¯}{X}}\right)\frac{\partial f_3}{\partial X_3}=Q\left({\displaystyle \underset{¯}{X}}\right)\frac{\partial g_3}{\partial X_3}{u}_{\sin }\Rightarrow u_{\sin }=\frac{\left(\frac{d}{dt}Q\left({\displaystyle \underset{¯}{X}}\right)-Q\left({\displaystyle \underset{¯}{X}}\right)\frac{\partial f_3}{\partial X_3}\right)}{\left(Q\left({\displaystyle \underset{¯}{X}}\right)\frac{\partial g_3}{\partial X_3}\right)}$ E82

The expression of u_sin determines how the dilution rate must be manipulated (varied) during the singular interval, being a function only of the state variables. When developed, the final expression obtained for u_sin is quite complex and extensive, which was implemented in computational programming language to evaluate the value of this variable over the singular interval.

The condition for stopping the integration of mass-balance equations during the period following the singular interval, which is conducted in batch mode (u = 0), is determined from the fact that for free final time (t_f) problems, how is the case, H(t_f) = 0. Thus

$H={{\displaystyle \underset{¯}{\lambda }}}^T\left({\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)+{\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)u\right){\displaystyle \underset{u=0}{\Rightarrow }}H={{\displaystyle \underset{¯}{\lambda }}}^T{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)\Rightarrow H={\lambda }_1{f}_1+{\lambda }_2{f}_2+{\lambda }_3{f}_3\Rightarrow H={\lambda }_1\frac{dX}{dt}+{\lambda }_2\frac{dS}{dt}+{\lambda }_3\frac{dP}{dt}$ E83

As the final conditions of the adjoint variables are λ₁ = 0, λ₂ = 0, and λ₃ = 1, the stopping condition at t = t_f is dP/dt = 0.

Thus, the problem-solving algorithm consisted of the following steps:

1. Integrate the mass-balance equations with u = 0 (batch operation) until the values of the state variables satisfy $Q\left({\displaystyle \underset{¯}{X}}\right)=0$ , being the instant that this occurs, the first switching time t₁;

2. From instant t₁, make u = u_sin until the reactor volume reaches its maximum value, which corresponds to the second switching time t₂;

3. Starting from time t₂, return the operation of the reactor with u = 0 until the stop condition is reached in t_f.

From the data reported by Costa [20] and summarized in Table 1, the mass-balance equations were numerically integrated to determinate t₁, i.e., the instant at which $Q\left({\displaystyle \underset{¯}{X}}\right)=0$ . According to the data presented in Table 2, this time instant can be determined as t₁ = 28.74 h, since $Q\left({\displaystyle \underset{¯}{X}}\right)$ changes its signal around this time value.

Using the computational program developed for the calculation of u_sin, it was verified that this control variable was practically constant during the singular interval, being 0.00213 h^-1 a value quite representative.

Kinetic parameters/operating variables	Values
μm (h^-1); B (g-S/g-X)	1.1×10^-1; 6×10^-3
πm (g g^-1 h^-1); km (g/L); ki (g/L); kh (h^-1)	4.0×10^-3; 1.0×10^-4; 1.0×10^-1; 1.0×10^-2
YX/S (g/g); YP/S (g/g); m(g g^-1 h^-1)	0.47; 1.2; 2.9×10^-2
X0 (g/L); S0 (g/L); P0 (g/L); Sf (g/L)	1.3; 69.0; 0.0; 500.0
V0 (L); Vf (L)	8.121; 10.0

Table 1.

Values of kinetic parameters and operating variables used in the case study of the optimization of penicillin production in fed-batch reactor (source: [20]).

With respect to the second switching time (t₂), this was determined from u_sin and the initial and final (maximum) volumes as follows:

$\frac{dV}{dt}=F=uV\Rightarrow \frac{dV}{V}=udt\Rightarrow \ln \left(\frac{V_f}{V_0}\right)=u\Delta t\Rightarrow \ln \left(\frac{10}{8.121}\right)=0.00213\Delta t\Rightarrow \Delta t=97.71\text{h}$ E84

Thus, the singular interval duration is of 97.71 h and the second switching time (t₂) is t₂ = t₁ + ∆t = 28.74 + 97.71 = 126.45 h. The mass-balance equations were then integrated until this time, and it was verified that, at the end of the integration, the stop condition (dP/dt=0) had been satisfied, dispensing with the complementary period of batch operation and reaching a final product concentration (P) of 6.35 g/L.

2.1.1. Simulation of the penicillin production bioprocess in fed-batch reactor under optimized conditions

The simulation of the bioprocess under optimized conditions was performed using a computer program in FORTRAN language. For the numerical integration of the ordinary differential equations corresponding to the mass balances, the variable-step fourth-order Runge-Kutta-Gill method was used [17]. The full profiles of the state variables during the bioprocess are shown in Figures 1 and 2.

Figure 1.

Temporal profiles of the S (dotted line), X (dashed line), and P (solid line) state variables during a fed-batch penicillin production bioprocess.

Figure 2.

Temporal profile of the state variable V (dotted line) and of the control variable u (solid line) during a fed-batch penicillin production bioprocess.

Due to the decoupling between biomass growth and product synthesis, this type of fermentation behaves as a biphasic process. Therefore, characteristic profiles of penicillin fermentation were obtained for the state variables as shown in Figure 1, i.e., a first phase of accumulation of the cell is observed in which the substrate is almost entirely consumed for this purpose, without associated product formation (trophophase). After this growth phase, the fed substrate is practically used for penicillin production since there is no further catabolic repression of the antibiotic synthesis due to the low levels of substrate concentration established in the reactor in this second phase (idiophase). In addition, the kinetic pattern observed is in agreement with that expected for a secondary metabolite, i.e., the production occurs mostly after cell growth.

Regarding the fermentation medium volume in the bioreactor, the behavior of this variable shown in Figure 2 was already expected since, during the batch operation, this volume is constant because there is no addition or removal of fermentation medium to or from the bioreactor. In the fed-batch operation with continuous feed of unfermented medium to the bioreactor at a constant flow rate, the volume increases linearly over time, as shown in Figure 2. The temporal profile exhibited by the control variable u (dilution rate) in a step format is due to the change in the bioreactor operation, from batch (u = 0) to fed-batch (u ≠ 0) mode.

2.2. Case study #2: determination of the optimal temperature profile in a batch bioreactor for penicillin production

In this case study, the fungal growth is described by the logistic law, a substrate-independent model for microorganism population dynamics. In addition, the production of penicillin is also modeled considering that the formation of antibiotic is not associated with cell growth, and that the product is degraded by hydrolysis according to a first-order kinetics. The mathematical model, comprising two ordinary differential equations corresponding to the mass balances of cell and product in a batch bioreactor, and containing four parameters, is represented by (more information about this model can be found in Ref. [21]):

$\frac{dX}{dt}=r_X={\mu }_{m}X\left(1-\frac{X}{X_m}\right)$ E85

$\frac{dP}{dt}=r_P-r_h=\beta X-k_{h}P$ E86

where t is the time, X is the cell concentration, P is the antibiotic concentration, r_X is the cell growth rate, r_P is the antibiotic production rate, r_h is the product’s hydrolysis rate, and μ_m, X_m, β, and k_h are the parameters of the model, according to the following meanings: μ_m is the maximum specific growth rate, X_m is the maximum possible cell concentration to be achieved, β is the constant of product formation not associated with growth, and k_h is the rate constant for the antibiotic hydrolysis reaction.

For the application of the Pontryagin’s maximum principle, the model variables were dimensionless and expressions describing the kinetic parameters (b_i) as a function of temperature (θ) were incorporated in order to extend the validity range of the model to non-isothermal conditions. These functions have shapes typical of those found in microbial or enzyme-catalyzed reactions (concave down parabolas). The hydrolysis of the antibiotic is neglected in this dimensionless version of the model, with this version given by the following equations [17]:

$\frac{d{y}_1}{d\tau }=b_1{y}_1-\frac{b_1}{b_2}{y}_1^2,y_1\left(0\right)=0.03$ E87

$\frac{d{y}_2}{d\tau }=b_3{y}_1,y_2\left(0\right)=0.0$ E88

where

• y₁ = dimensionless concentration of cell (-); y₂ = dimensionless concentration of product (-); τ = dimensionless time, 0 ≤ τ ≤ 1 (-)

  $b_1=w_1\text{[}\frac{1.0-w_2{\left(\theta -w_3\right)}^2}{1.0-w_2{\left(25-w_3\right)}^2}\text{]};b_2=w_4\text{[}\frac{1.0-w_2{\left(\theta -w_3\right)}^2}{1.0-w_2{\left(25-w_3\right)}^2}\text{]};b_3=w_5\text{[}\frac{1.0-w_2{\left(\theta -w_6\right)}^2}{1.0-w_2{\left(25-w_6\right)}^2}\text{]}$ E89-91

• $w_1=13.1;w_2=0.005;w_3=30°\text{C}$

• $w_4=0.94;w_5=1.71;w_6=20°\text{C}$

As in this case, ${\displaystyle \underset{¯}{g}}\left({\displaystyle \underset{¯}{X}}\right)=0$ and u = 0, due to the reactor being operated in a batch mode, the mass-balances equations are simplified to

$\frac{d{\displaystyle \underset{¯}{X}}}{dt}={\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)$ E92

where

${\displaystyle \underset{¯}{X}}=\left[\begin{array}{c}{y}_1\\ y_2\end{array}\right];{\displaystyle \underset{¯}{f}}\left({\displaystyle \underset{¯}{X}}\right)=\left[\begin{array}{c}{f}_1\\ f_2\end{array}\right]=\text{[}\begin{array}{c}{b}_1{y}_1-\frac{b_1}{b_2}{y}_1^2\\ b_3{y}_1\end{array}\text{]}$ E93-94

As previously established in the first case study, the Hamiltonian is given by

$\begin{array}{l}H={\underset{¯}{\lambda }}^T\left(\underset{¯}{f}\left(\underset{¯}{X}\right)+\underset{0}{\underbrace{\underset{¯}{g}\left(\underset{¯}{X}\right)u}}\right)\Rightarrow H={\underset{¯}{\lambda }}^T\underset{¯}{f}\left(\underset{¯}{X}\right);{\underset{¯}{\lambda }}^T=\text{[}{\lambda }_1{\lambda }_2\text{]}\Rightarrow \\ H=\text{[}{\lambda }_1{\lambda }_2\text{]}\left[\begin{array}{l}{f}_1\\ f_2\end{array}\right]={\lambda }_1{f}_1+{\lambda }_2{f}_2={\lambda }_1\left(b_1{y}_1-\frac{b_1}{b_2}{y}_1^2\right)+{\lambda }_2\left(b_3{y}_1\right)\end{array}$ E95

The temporal variation rates of the adjoint variables λ₁ and λ₂ are formulated as

$\frac{d{\displaystyle \underset{¯}{\lambda }}}{d\tau }=-\frac{\partial H}{\partial {\displaystyle \underset{¯}{X}}}\Rightarrow \frac{d}{d\tau }\left[\begin{array}{l}{\lambda }_1\\ {\lambda }_2\end{array}\right]=\text{[}\begin{array}{l}-{\lambda }_1{b}_1+2\frac{b_1}{b_2}{y}_1-{\lambda }_2{b}_3\\ 0\end{array}\text{]}$ E96

From the previous equation, the following equations can be derived

$\frac{d{\lambda }_1}{d\tau }=-b_1{\lambda }_1+2\frac{b_1}{b_2}{y}_1{\lambda }_1-b_3{\lambda }_2$ E97

$\frac{d{\lambda }_2}{d\tau }=0$ E98

The necessary condition for the optimization of the bioprocess is

$\frac{\partial H}{\partial \theta }=0\Rightarrow \frac{\partial H}{\partial \theta }={\lambda }_1\left(y_1\left(\frac{\partial b{}_1}{\partial \theta }\right)-y_1^2\frac{\partial \left(b_1/b_2\right)}{\partial \theta }\right)+{\lambda }_2\left(y_1\frac{\partial b_3}{\partial \theta }\right)=0$ E99

From the expressions b_i = b_i(θ), the following derivatives can be obtained

$\frac{\partial b{}_1}{\partial \theta }=-\text{[}\frac{2w_1{w}_2\left(\theta -w_3\right)}{1.0-w_2{\left(25-w_3\right)}^2}\text{]}\text{;}\frac{\partial \left(b_1/b_2\right)}{\partial \theta }=0;\frac{\partial b_3}{\partial \theta }=-\text{[}\frac{2w_5{w}_2\left(\theta -w_6\right)}{1.0-w_2{\left(25-w_6\right)}^2}\text{]}$ E100-102

By inserting the derivatives of the parameters with respect to the temperature into the expression of ∂H/∂θ = 0, the expression of the optimal temperature profile (θ_opt) is obtained as follows:

${\theta }_{opt}=\text{[}\frac{2{\lambda }_1{y}_1{w}_1{w}_2{w}_3}{1.0-w_2{\left(25-w_3\right)}^2}+\frac{2y_1{w}_5{w}_2{w}_6}{1.0-w_2{\left(25-w_6\right)}^2}\text{]}/\text{[}\frac{2{\lambda }_1{y}_1{w}_1{w}_2}{1.0-w_2{\left(25-w_3\right)}^2}+\frac{2y_1{w}_5{w}_2}{1.0-w_2{\left(25-w_6\right)}^2}\text{]}$ E103

As previously demonstrated, when the objective is to maximize the antibiotic concentration at the end of the bioprocess, it is necessary that λ₁(1) = 0 and λ₂(1) = 1. Since dλ₂/dt = 0, the second condition requires that λ₂ be constant and equal to 1.0 in the entire time domain, i.e., λ₂ = 1.0 for 0 ≤ τ ≤ 1.

Several numerical methods have been developed to solve this two-point boundary-value problem arising from the application of the maximum principle of Pontryagin to a batch penicillin production bioprocess. Constantinides and Mostoufi [17] used the orthogonal collocation method to solve this problem, justifying that this method is more accurate than the finite difference method. The problem was solved here using a much simpler numerical method than that of the orthogonal collocation to integrate the differential equations, which is the variable-step fourth-order Runge-Kutta-Gill method [17]. Thus, the algorithm for solving the problem consisted of the following steps:

1. Assignment of an initial value for λ₁(0);

2. Integration of the system of ODEs from τ = 0 to τ = 1 and verification if λ₁(1) = 0. If not, assign a new value to λ₁(0) until the final condition is satisfied.

In order to make the computational algorithm autonomous for the determination of λ₁(0), the Newton-Raphson method [17] was coupled to the numerical integration method, solving the following non-linear algebraic equation:

$r\left({\lambda }_1\left(0\right)\right)={\text{[}{\lambda }_1\left(1\right)\text{]}}_{calculated}-{\displaystyle \underset{0}{{\displaystyle \underbrace{{\text{[}{\lambda }_1\left(1\right)\text{]}}_{specified}}}}}=0\Rightarrow r\left({\lambda }_1\left(0\right)\right)={\text{[}{\lambda }_1\left(1\right)\text{]}}_{calculated}=0$ E104

2.2.1. Simulation of the penicillin production bioprocess in batch reactor under optimized non-isothermal conditions

The proposed algorithm was implemented in FORTRAN programming language, and the profiles of the state variables (y₁, y₂, and θ) are presented in Figures 3–5. These profiles are in strict agreement with those reported by Constantinides and Mostoufi [17] when they used the orthogonal collocation method to solve this problem. The value determined for λ₁ at τ = 0 was λ₁(0) = 3.61210035.

Figure 3.

Cell dimensionless concentration profile during a non-isothermal penicillin fermentation.

Figure 4.

Product dimensionless concentration profile during a non-isothermal penicillin fermentation.

Figure 5.

Exact and approximate optimal temperature profiles for a non-isothermal penicillin fermentation.

The cell concentration profile shown in Figure 3 depicts the main phases involved in a typical microbial growth curve, i.e., the exponential, stationary, and decline phases. The decline phase is attributed to the negative effects of low temperatures on cell growth. In Figure 4, concerning the penicillin production dynamics, an initial short lag phase can be observed, followed by a transition phase in which penicillin production is initiated, until a final linear production phase is achieved.

According to the presented formulation (bioprocess model and Pontryagin’s maximum principle), the optimum temperature profile varies between 20 and 30°C following the curve (a) shown in Figure 5. This profile suggests a variable operating temperature during the growth and penicillin production phases, contradicting the standard industrial operating procedure of maintaining a constant temperature throughout the bioprocess. Particularly during the penicillin production phase, the profile prescribes a decrease in the operating temperature so that a high concentration of antibiotic is reached at the end of the bioprocess. The temperature profile (a) shown in Figure 5 may bring some practical difficulty to its programming/execution. In this context, an approximate profile derived from the exact profile, such as that represented by the curve (b) in Figure 5, may make the temperature programming strategy more feasible. This proposal is in agreement with that reported by Bailey and Ollis [22], i.e., the temperature schedule predicted by these calculations can be closely approximated in industrial practice with little added cost. The approximate profile was built from the following equations, which were based on the analysis of the temperature data generated by the exact profile:

• $0.0\le \tau \le 0.4\text{: }\theta (°\text{C})={\theta }_{(\tau =0.0)}+\left(\frac{{\theta }_{(\tau =0.4)}-{\theta }_{(\tau =0.0)}}{0.4-0.0}\right)\left(\tau -0.0\right)\Rightarrow \theta =29.65-11.625\tau$

• $0.4<\tau <0.9\text{: }\theta (°\text{C})=25.0$

• $0.9\le \tau \le 1.0\text{: }\theta (°\text{C})={\theta }_{(\tau =0.9)}+\left(\frac{{\theta }_{(\tau =1.0)}-{\theta }_{(\tau =0.9)}}{1.0-0.9}\right)\left(\tau -0.9\right)\Rightarrow \theta =25.0-50\left(\tau -0.9\right)$

2.2.2. Simulation of the penicillin production bioprocess in batch reactor under non-optimized isothermal conditions

A pertinent simulation to be performed is one under isothermal conditions to verify whether this thermal operation mode is, in fact, less productive in penicillin than non-isothermal mode following an optimal temperature profile. For this purpose, simulations were performed for constant temperatures of 20, 25, and 30°C and the results were compared with those obtained with the optimized temperature profile (Figure 6). Figure 6(a) illustrates the well-known fact that high temperatures (30°C) favor the growth of the fungus, while low temperatures (20°C) favor the synthesis of the antibiotic since relative to the amount of penicillin produced, more and less biomass was accumulated at these respective temperature levels [22]. It is observed in Figure 6(b) that the isothermal operation at an intermediate temperature to those investigated (θ = 25°C) presents a performance very similar to the operation with optimized temperature profile, becoming a viable alternative if the variable temperature strategy cannot be implemented, although with lower productivity, according to the data presented in Table 3. Despite the fact that operation at a fixed temperature between 24 and 25°C predominates in current industrial practice, the benefits of temperature programming during batch antibiotic fermentations are clear.

t (h)	X (g/L)	S (g/L)	P (g/L)	Q(X)
28.742	30.09	4.16×10^-3	1.39×10^-2	125799.24
28.743	30.09	3.06 ×10^-3	1.40×10^-2	−86336.84

Table 2.

Data of the numerical integration of the mass-balance equations used to determine the first switching time (t₁) during the penicillin production in fed-batch reactor.

θ (°C)	y1 (τ =1)	y2 (τ =1)	(y2/y1)τ =1
θopt	0.82	1.22	1.49
20	0.53	0.65	1.23
25	0.94	1.18	1.25
30	1.07	0.80	0.75

Table 3.

Final values of the dimensionless concentration of cells (y₁) and product (y₂) in isothermal and non-isothermal penicillin fermentations.

Figure 6.

Dimensionless concentration profiles of cell and product during isothermal penicillin fermentations at different temperatures.

3. Conclusions

In this chapter, the usefulness of the Pontryagin’s maximum principle has been demonstrated for the optimization and operation of complex antibiotic production bioprocesses such as those conducted in batch and fed-batch reactors under isothermal/non-isothermal conditions. By applying this principle, it was possible to determine the optimal profile of temperature in batch reactors and substrate feed rate in fed-batch reactors that maximize the antibiotic concentration at the end of the bioprocess. Although having a rather complex mathematical formulation, the Pontryagin’s maximum principle can be classified as a powerful and suitable tool for the optimization, control, and model-driven operation of bioprocesses aiming at maximum productivity of bioproducts. However, for the application of this principle, it is necessary to dispose a mathematical model, preferably phenomenological and representative of the bioprocess, in order to evaluate whether or not the solution found for a given problem is feasible. In the present study, two classical phenomenological models of penicillin production bioprocesses were used, together with the Pontryagin’s maximum principle, aiming to determine the optimal operating conditions for the production of antibiotic, and the solutions found are considered feasible and can be implemented in real cases. However, a more complete mathematical model, incorporating the medium oxygenation state, could provide better bioprocess control, since the productivity in penicillin fermentations is highly dependent upon dissolved oxygen concentration, with its critical level being around 30% of saturation. In the models used here, the dissolved oxygen concentration was implicitly assumed to be non-limiting of the bioprocess, making this a rather restrictive hypothesis.

Acknowledgments

The author wishes to thank CNPq for their financial support (protocol number: 455487/2014-6).