Dynamic Model for Biomass and Proteins Production by Three Bacillus Thuringiensis ssp Kurstaki Strains

: Bacillus thuringiensis is a microorganism used for the production of biopesticides worldwide. In the present paper, different kinetic models were analyzed to study and compare three different strains of Bt ssp kurstaki (LIP, BLB1, and HD1). Bioperformances (vegetative cell, spore, substrate, and protein) and successive culture phases (oxidative growth, limitation and sporulation, and protein release) were depicted with an overarching aim to estimate total protein productivity, yield, and titer. In the end, two models were calibrated using experimental dataset (11 batches culture in 3 L bioreactor with semisynthetic medium), subsequently validated, and statistically compared. Both models satisfactorily followed the dynamics of the experimental data. Finally, a dynamic model was selected following the Akaike information criterion (AIC).


Introduction
B. thuringiensis is a facultative anaerobic Gram-positive sporulating bacterium, frequently used in the production of some biopesticides and as a source of genes for transgenic expression in plants [1]. It usually inhabits different environments, among which soil, settled dust, insects, water, and others have been identified [2]. B. thuringiensis has been shown to be toxic to various organisms, such as lepidopterans, coleopterans, dipterans, or nematodes, but is considered safe for mammals. Thus, the products based on B. thuringiensis (Bt) provide effective and environmentally benign control of several insects in agricultural, forestry, and disease-vector applications [3]. This insecticidal activity is mainly due to the production of some intracellular inclusions (called σ-endotoxins) during the sporulation phase of B. thuringiensis cells.
Most of the biopesticides distributed in the world are mainly based on Btk. HD1 strain. However, two recent strains, identified as Btk. LIP (from Lebanese soil), and BLB1 (from Tunisian soil), have been isolated and described to be more efficient than HD1 [4], and, therefore, will be studied in this work.
Due to the several changes of cell physiology during the σ-endotoxins production bioprocess (exponential growth, formation of inclusion, formation of spore, lysis), B. thuringiensis culture is considered a laborious process. Although one possibility to optimize B. thuringiensis culture is through mathematical models, there are not many mathematical models that describe the dynamics of the growth phases of B. thuringiensis culture [5]. Holmberg and Sievanen [6] proposed a model based on Monod kinetics to describe the 2 of 16 relationship between cell growth and toxin production. Later, Rivera et al. (1999) [7] showed the Monod model limitations when they tried to describe the biomass diversity of B. thuringiensis. In their model, they assumed the presence of two kind of cells in the B. thuringiensis culture; those available to multiply, and those that have become spores. Thus, they divided them between the biomass of the vegetative cells and the biomass of the sporeforming cells. In addition, they used the Monod model to describe the relationship between vegetative cell growth and substrate concentration, as did Holmberg and Sievanen [6].
Furthermore, Popovic et al. [8] proposed a model that considered a minimum level of poly-β-hydroxy butyric acid (PHB) required in cells at the beginning of sporulation for efficient sporulation and endotoxins productions. Additionally, they used Contois kinetics to describe the growth of the cells, considering that this model fits better to the experimental data than the expression of Monod.
Therefore, this work proposes a dynamic model for B. thuringiensis. Section 2 describes experimental data and the dynamic model for B. thuringiensis. The simulation results and the performance evaluation are shown in Section 3. Finally, conclusions and perspectives close this paper.

Materials and Methods
The materials and methods used to generate the set of experimental data are described in this section. Then, experimental data are introduced, and assumption and formulation of models are explained.

Microorganism and Culture Media
Three B. thuringiensis ssp. kurstaki strains were used in the present work: a Lebanese strain LIP [4], a Tunisian strain BLB1 [9], and HD1 strain used as a reference (industrial gold standard) [10]. Luria broth (LB) medium was used for inoculum production, whereas a semisynthetic medium (SSM), defined by Sarrafzadeh et al. [11], was used for fermentation assays. Their compositions (g·L -1 ) are described in Table 1. For the SSM, concentrated glucose (Sol 2) and all salts solutions (Sol 3-5) were prepared and sterilized separately and added before inoculation to the rest of the medium (Sol 1) previously sterilized. corresponding to an initial OD600 = 0.15 were used to inoculate 1 L Erlenmeyer flasks containing 100 mL LB medium. After 10-12 h of incubation at 30 • C, in a rotary shaker set at 200-230 rpm, the OD600 was determined. The culture broth was used to inoculate the bioreactor containing the studied media to start with an initial OD600 of 0.15.

Culture Conditions
Several fermentations were conducted over 48 h in batch mode at 30 • C in a 3 L Biostat B plus fermenter (Sartorius, Göttingen, Germany) containing 1.8 L of the SSM medium and with continuous regulation of pH at 6.8 using 1 M H 2 SO 4 and 3 M NaOH. Dissolved oxygen was continuously monitored by an optical oxygen sensor and maintained at 25% pO2-sat with a constant aeration rate (VVM = 10 with Qair = 0.18 min·L −1 ) and variable stirring (from 250 to 1200 rpm). Foaming was controlled using an antifoam (Emultrol DFM DV-14 FG), through the fermentation process.

Analysis and Sampling Strategy
Several samples were collected from the Bt broth during experiments, and substrate, biomass, and product analyses (glucose, cell and spore counting, protein) were conducted to determine biokinetics.

Detection of Sporulation
The diverse cell states were distinguished, during the fermentation process, based on their morphological differences and the refractile nature of the endospores, using a phase contrast microscope (ZeissPrima Pro, Paris, France, ×100 oil).

Biomass Analyses Cells and Spores Counts
The follow-up of the biopesticides production was checked by estimating viable cell counts (VC) and spore counts (SC) by plate counts. To determine VC and SC, the withdrawn samples were serially diluted, spread on LB plates and incubated at 30 • C for 16-18 h. For determining SC, the appropriately diluted samples were heated at 85 • C for 15 min and cooled for 5 min before spreading onto LB plates. Number of colonies should be between 20 and 300 to be acceptable. All analyses were realized in triplicate.

Cell Dry Weight
A known amount of sample (1 to 20 mL) was filtered via nitrocellulose membrane (0.2 µm) and the membrane was then dried at 70 • C (24 h). Biomass dry weight is determined by differential weighing of the filter before the filtration and after filtration and drying.

Quantification of Proteins Production
In order to estimate the concentration of total proteins (mainly composed of δendotoxin) produced during the fermentation, 1 mL sample was centrifuged at 13,000 rpm for 5 min at 4 • C. The supernatant was collected for other analysis and the pellet was washed twice with cold NaCl 1M and four times with cold water. The protein crystal was then dissolved in 0.05 N NaOH for 2 to 3 h at 30 • C in a rotary shaker (200 rpm). The suspension was then centrifuged at 13,000 rpm for 5 min, and the supernatant containing the solubilized proteins was recuperated.
The concentration of the proteins in this supernatant was determined by Bradford assay [12] using bovine serum albumin (BSA) as a protein standard. Absorbances were measured after 10 min at 595 nm (2300 EnSpire Multilabel Plate Reader). The obtained value was the average of three measures of the same sample (microwell plate). Considering our protocol, protein concentration estimates the total protein production after separation but not specifically the δendotoxin, even if it is the dominant fraction.

Sugar Analysis
The sugar concentrations were determined using HPLC-UV. The HPLC assays were performed using an Ultimate 3000 RSLC/MWD/RI/CAD. A mobile phase of 5 mM H 2 SO 4 with a flow rate of 0.6 mL·min −1 was used. The mobile phase was filtered and degassed through a 0.2 µm cellulose nitrate membrane. The samples and standards were also filtered before injection into the HPLC.

Experimental Data
Between three to four batch cultures per strain were carried out ( Table 2). Two batches per strain were used to perform parameter calibration, and between one and two batches were used to validate the models. The experimental datasets for each strain are presented in Figures 1-3. It is relevant to indicate that in batch 07, dry matter measurement was estimated by OD600nm for exponential growth phase.

Btk. HD1 Strain
Similar data and milestones were obtained for the HD1 strain ( Figure 2). Around 10 h, maximum biomass concentration and substrate depletion were reached. However, in batch 07, no values were recorded for the biomass concentration after exponential growth phase due to technical misplaced measurements. After 20 h, sporulation rate was achieved, and protein content plateaued after 30 h. Final protein content was around 0.8-1 g/L, except for batch 07 (0.4 g/L).   Figure 3 presents identical variables and leads to identify the same milestones and critical time as described above with Btk BLB1 and HD1 strains. LIP strain exhibited the lowest protein concentration, close to 0.2 g/L. This result could be explained by a lower cell lysis rate; therefore, protein crystals were not released in the supernatant.

Model Assumptions
The main features of dynamic models include the key parameters describing bioperformances (vegetative cell, spores, substrate, proteins) and associated kinetics, consider- During all cultures, several milestones should be identified: (i) the maximum cell biomass production, (ii) substrate depletion, (iii) full sporulation informing about proteins production, and (iv) its release in supernatant due to full cell lysis.
For the BLB1 strain, data were taken from four batches (01, 02, 05, and 06). Figure 1 presents the graphical experimental data. Batch 11 was used to validate the model.
The maximum biomass concentration was reached after approximately 10 h culture for the four batches, and then began to decrease. As expected, this time approximately coincided with the substrate depletion, corresponding to a glucose concentration close to 0 g·L −1 . In addition, the concentration of spores began to increase at this moment, since the limitation of the substrate induced their formation. After 20 h, sporulation reached a plateau value. Cell lysis was fully achieved after 30 h, as indicated by protein release into supernatant. Finally, protein content reached around 0.8-1 g·L −1 with BLB1 stain.

Btk. HD1 Strain
Similar data and milestones were obtained for the HD1 strain ( Figure 2). Around 10 h, maximum biomass concentration and substrate depletion were reached. However, in batch 07, no values were recorded for the biomass concentration after exponential growth phase due to technical misplaced measurements. After 20 h, sporulation rate was achieved, and protein content plateaued after 30 h. Final protein content was around 0.8-1 g/L, except for batch 07 (0.4 g/L). Figure 3 presents identical variables and leads to identify the same milestones and critical time as described above with Btk. BLB1 and HD1 strains. LIP strain exhibited the lowest protein concentration, close to 0.2 g/L. This result could be explained by a lower cell lysis rate; therefore, protein crystals were not released in the supernatant.

Model Assumptions
The main features of dynamic models include the key parameters describing bioperformances (vegetative cell, spores, substrate, proteins) and associated kinetics, considering successive phases (oxidative growth, limitation, sporulation, and protein release), during bioproduction. The mass balance equations on each compound are shown in Equations (1), (2), (4)-(7). Equation (1) represents the evolution of biomass concentration with respect to time, while the relationship between bacterial growth and substrate consumption is shown in Equations (1) and (2).
where X is the biomass concentration (g·L −1 ), S is the concentration of substrate (g·L −1 ), Y1 is the yield coefficient between the biomass and the substrate (gBiomass/gGlucose), and k d is the death rate (h −1 ). The cell growth process is represented by the Contois expression, as follows: where µ is the specific growth rate (h −1 ) and µ max is the maximum specific growth rate (h −1 ), a constant defined for a substrate concentration; X1 is the concentration of biomass (g·L −1 ); S1 is the concentration of glucose (g·L −1 ); and Kc is a saturation constant. Equations (4)- (7) show the mass balance for proteins and spores. In the first model, proteins and spores are correlated with biomass (Model 1). In the second model, α and β parameters were used to unassociated proteins and spores from biomass production (Model 2). The corresponding models are shown in Equations (4) and (5) (Model 1) and Equations (6) and (7) (Model 2).
SSR is sum of squared regression, and SST is sum of squared total.
where n represents the number of data, p the number of parameters, W the weighting matrix, and X and X are the data and estimated data, respectively [13].
The main parameter used to determine the model that best fits the data is the AICc parameter [13]. The AIC criterion is one of the most popular for the comparison of models because it considers the number of parameters, the number of data, and the residuals, making it a parameter that balances the complexity of the model and the fit of the data [14]. Additionally, the parameter correction (AICc) gives accurate results for a larger number of datasets. Thus, the model with the lowest value for AICc is selected to represent the experimental data more adequately.
Parameter calibration was carried out in MATLAB using a particle swarm optimization (PSO) algorithm. This method, as its name implies, is inspired by the behavior of swarms of insects in nature. Thus, for a set of variables to be optimized, the method begins by placing random particles in the search space, but then a series of rules are established considering each parameter and the set of parameters ("swarm") globally. Thus, the variables are optimized quite well, and few computational resources are spent, becoming a fast method in convergence, and simple in application [15].

Results
This section presents the results obtained through various simulations carried out in MATLAB. It is divided into two main sections: model calibration and model validation. Initially, the parameters of each model were estimated (six for Model 1 and eight for Model 2) with the experimental data of two batches for each strain. Subsequently, the parameters found were used to observe the behavior of the four state variables (biomass, glucose, proteins, and spores concentrations) in two batches for the BLB1 strain and one batch for the HD1 and LIP strains. Likewise, the results obtained were compared with the experimental data. In this way, the validation of the parameters found in the calibration phase was carried out. Moreover, to compare models 1 and 2, a series of statistical parameters were calculated from which the selection of the model that best fits the experimental data is facilitated.

Model Calibration BLB1 Strain
Kinetics parameters of the B. thuringiensis culture were calibrated for three strains: BLB1 (Table 3), HD1 (Table 5), and LIP ( Table 7). The maximum specific growth rate (µmax) was between 1.15 h−1 for the BLB1 strain (Model 1) and 0.39 h −1 for the LIP strain (Model 2). These results are within ranges similar to those reported by Holmberg and Sievanen (1980), who reported values between 1.90 and 0.17 h −1 [6], and Atehortúa et al. (2007), who reported values between 0.80 and 0.58 h −1 [16]. The death rate (k d ) was between 0.0458 (BLB1 strain) and 0.0184 h −1 (LIP strain), which coincided with previous results in the literature (between 0 and 0.13 h −1 ) [6]. The yield coefficient between the biomass and the substrate (Y1) was between 0.49 (gBiomass·gGlucose −1 ) and 0.96 (gBiomass·gGlucose −1 ) for all strains. The comparison between Models 1 and 2 and the experimental data for the BLB1 strain are shown in Figure 4. According to the figure, Models 1 and 2 did not have very noticeable differences; therefore, the alpha and beta parameters of Model 2 did not have a great impact on the modeling. Both models showed a satisfactory fit to the experimental data; however, the statistical study will give precise information on the best model. It is important to note that the quantification of the spore concentration and protein concentration are more subject to systematic error than biomass and glucose, which may explain the discrepancies between models and measurements.
Model 1 had higher values for µmax and Kc, although other parameters remained with similar values. Furthermore, for the constants alpha and beta, very low values of 0.0001 were obtained, which confirmed that both models are similar. Table 4 presents the statistical coefficients that allowed comparing Model 1 and Model 2 for BLB1 strain. The statistical coefficients showed a good fit to the experimental data. Spores concentration was the variable with the worst fit. These results are supported by Figure 4, since the furthest experimental values of the two models can be observed in it. strain are shown in Figure 4. According to the figure, Models 1 and 2 did not have very noticeable differences; therefore, the alpha and beta parameters of Model 2 did not have a great impact on the modeling. Both models showed a satisfactory fit to the experimental data; however, the statistical study will give precise information on the best model. It is important to note that the quantification of the spore concentration and protein concentration are more subject to systematic error than biomass and glucose, which may explain the discrepancies between models and measurements.   As far as the determination coefficient (R2), values greater than 0.98 were observed for glucose measurements, which indicates that this variable has the best fit. However, as mentioned above, the parameter of greatest interest is the AICc since it takes into account several important aspects. The model that presented the lowest AICc values for the four variables was Model 1. This model does not include any extra constants, which makes it a less complex model than Model 2, but also predicts the behavior of the variables studied.

Model Calibration HD1 Strain
The results for the HD1 strain are shown in Figure 5 and Table 5. It was shown that there are big differences between Model 1 and Model 2. Protein concentration and spore concentration were the variables in which these differences were most visible according to Figure 5, which makes sense since the alpha and beta constants present in Model 2 have a direct influence on these two variables. Additionally, both models presented a very good fit for the biomass and substrate concentration. Processes 2021, 9, x FOR PEER REVIEW 11 of 18 The optimized parameters showed higher values in Model 1 than in Model 2 for almost all parameters. Although the values of µ max were lower than those found for the BLB1 strain, parameters such as Y1 and alpha and beta constants were higher than those obtained with the BLB1 strain. The results of the statistical parameters for the calibration of the HD1 strain are presented in Table 6. In a similar way to the BLB1 strain, the variable that best fits the models according to the coefficient of determination was glucose, even reaching a value of 1 for batch 4 and Model 2. In general, the R2 and RMSE coefficients showed a better fit of the models with the experimental data for the HD1 strain than the previously analyzed BLB1 strain. In fact, for HD1 strain, the experimental data of the spore concentration have a better fit than those obtained for the BLB1 strain.
Since for the HD1 strain, Model 1 showed the lowest AICc values, this model was considered as the most appropriate to predict the behavior of the HD1 strain according to the results of the calibration. This means that Model 1 showed the best parsimony.  The optimized parameters showed higher values in Model 1 than in Model 2 for almost all parameters. Although the values of µmax were lower than those found for the BLB1 strain, parameters such as Y1 and alpha and beta constants were higher than those obtained with the BLB1 strain.
The results of the statistical parameters for the calibration of the HD1 strain are presented in Table 6. In a similar way to the BLB1 strain, the variable that best fits the models according to the coefficient of determination was glucose, even reaching a value of 1 for batch 4 and Model 2. In general, the R2 and RMSE coefficients showed a better fit of the models with the experimental data for the HD1 strain than the previously analyzed BLB1 strain. In fact, for HD1 strain, the experimental data of the spore concentration have a better fit than those obtained for the BLB1 strain.
Since for the HD1 strain, Model 1 showed the lowest AICc values, this model was considered as the most appropriate to predict the behavior of the HD1 strain according to the results of the calibration. This means that Model 1 showed the best parsimony. Figure 6 shows the results obtained for the LIP strain. Graphically, Model 1 and Model 2 showed great similarities except for protein concentration. The values obtained for the parameters and the statistical coefficients reflect these differences.   Figure 6 shows the results obtained for the LIP strain. Graphically, Model 1 and Model 2 showed great similarities except for protein concentration. The values obtained for the parameters and the statistical coefficients reflect these differences.  Table 7 summarizes the parameter values of both models. It is noteworthy that the proteins/substrate yield coefficient (Y2) showed very low values.

Model Calibration LIP Strain
Moreover, the calibrated parameters showed a higher value for the alpha parameter than the beta one. Therefore, the beta parameter, has a very small value and little influence on Model 2.   Table 7 summarizes the parameter values of both models. It is noteworthy that the proteins/substrate yield coefficient (Y2) showed very low values. Moreover, the calibrated parameters showed a higher value for the alpha parameter than the beta one. Therefore, the beta parameter, has a very small value and little influence on Model 2. Table 8 indicates the values of the statistical parameters for the LIP strain. The statistical parameters showed a good fit of the models, presenting very low values of the determination coefficient only for the protein concentration in batch 10. Glucose continues to be the variable that has the best fits between two models. Additionally, as in the BLB1 and HD1 cases, Model 1 obtained the lowest values for AICc, making it the model that could best predict the behavior of the LIP strain.

Model Validation
Several datasets of each strain, different from those used in the calibration of the models, were used to validate the results obtained previously. Batch 2 and 5 were used to validate the parameters obtained for the BLB1 strain, the data from batch 7 were used for the HD1 strain, and, finally, batch 8 helped validate the parameters of the LIP strain. Figure 7 shows the results of the validation for the BLB1 strain, and Table 9 shows the respective statistical coefficients. For both experiments, the calibrated parameters fit the experimental data very well. According to Figure 7, Models 1 and 2 behaved similarly and there were no noticeable differences.  Statistical analysis showed that Model 1 fit better than Model 2 because it has the lowest values of the AICc coefficient. As seen graphically, the statistical parameters showed less adjustment for some variables of batch 02 than for batch 05. Data from batch 7 were used for validation of the parameters obtained in HD1 strain calibration. The results are shown in Figure 8 and Table 10. As said before, no values were recorded for the biomass concentration after exponential growth phase due to technical misplaced measurements, which is reflected in Figure 8. However, simulations showed that biomass during exponential growth phase and the other state variables fit adequately. The set of statistical coefficients that express the effectiveness of the models is expressed in Table 10. Statistical analysis showed that Model 1 fit better than Model 2 because it has the lowest values of the AICc coefficient. As seen graphically, the statistical parameters showed less adjustment for some variables of batch 02 than for batch 05.

HD1 Strain
Data from batch 7 were used for validation of the parameters obtained in HD1 strain calibration. The results are shown in Figure 8 and Table 10. As said before, no values were recorded for the biomass concentration after exponential growth phase due to technical misplaced measurements, which is reflected in Figure 8. However, simulations showed that biomass during exponential growth phase and the other state variables fit adequately. The set of statistical coefficients that express the effectiveness of the models is expressed in Table 10. Both models fit quite well for glucose concentration and followed the dynamics of spores concentration. The statistical coefficients showed the worst results for spores. As demonstrated previously, according to the AICc criterion, Model 1 should be the one used to represent the data of the HD1 strain.  Both models fit quite well for glucose concentration and followed the dynamics of spores concentration. The statistical coefficients showed the worst results for spores. As demonstrated previously, according to the AICc criterion, Model 1 should be the one used to represent the data of the HD1 strain. Figure 9 show the results of the validation for the LIP strain. Similar to batch 7 (HD1 strain), batch 8, which corresponds to the LIP strain, showed a biomass measurement problem. However, glucose concentration was well represented by the two models. Although the models followed the dynamics for the concentration in proteins and spores, a slight lag was evident for the spores.  Table 11 presents results for statistics parameters in LIP validation. As for the two strains, the glucose data showed the best fit. However, the coefficient of determ tion for proteins and spores showed a good fit of the models. As in all the cases stu in this report, the most suitable model to predict and represent the experimental d Model 1, according to the AICc criterion.   Table 11 presents results for statistics parameters in LIP validation. As for the other two strains, the glucose data showed the best fit. However, the coefficient of determination for proteins and spores showed a good fit of the models. As in all the cases studied in this report, the most suitable model to predict and represent the experimental data is Model 1, according to the AICc criterion.

Conclusions
The objectives of approach were to model bioperformances (vegetative cell, spore, substrate, and protein) considering different B. thuringiensis ssp. kurstaki strains and successive culture phases (oxidative growth, limitation and sporulation, protein release). Our overarching aim to estimate total proteins production (mainly composed of δ-endotoxin) was successfully achieved. Initially, the bibliographic research allowed understanding of the context and the different phenomena involved in the study of the B. thuringiensis culture, such as the particular life cycle of these microorganisms and the importance of the endotoxins produced.
B. thuringiensis is an important microorganism for the biopesticide market worldwide. The experimental simulations developed in the present study and based on B. thuringiensis cultures made it possible to analyze the behavior of the concentration in biomass, substrate, proteins, and spores and adjust two models to the experimental datasets. The calibration of both models allowed to calculate the kinetic parameters of the culture, and the experimental data presented a good fit. Likewise, the models were validated in a satisfactory way.
For the selection of the best model, the AICc criterion was used, which, for all batches, showed better results for Model 1 due to its parsimony. Additionally, although the BLB1 strain showed the highest maximum specific growth rate (µmax), the HD1 strain presented the highest biomass/substrate yield coefficient values (Y1), in opposition to the LIP strain which presented the lowest values for this yield. As for the production of proteins, mainly used for insecticidal toxicity, the BLB1 strain presented the highest concentration and proteins/substrate yield coefficient (Y2), while the LIP strain showed the lowest values for this yield.