*Article* **Advanced Kinetic Modeling of Bio-co-polymer Poly(3-hydroxybutyrate-***co***-3-hydroxyvalerate) Production Using Fructose and Propionate as Carbon Sources**

**Stefanie Duvigneau 1,\*, Robert Dürr 2, Jessica Behrens 1,2 and Achim Kienle 1,2**


**Abstract:** Biopolymers are a promising alternative to petroleum-based plastic raw materials. They are bio-based, non-toxic and degradable under environmental conditions. In addition to the homopolymer poly(3-hydroxybutyrate) (PHB), there are a number of co-polymers that have a broad range of applications and are easier to process in comparison to PHB. The most prominent representative from this group of bio-copolymers is poly(3-hydroxybutyrate-*co*-3-hydroxyvalerate) (PHBV). In this article, we show a new kinetic model that describes the PHBV production from fructose and propionic acid in *Cupriavidus necator* (*C. necator*). The developed model is used to analyze the effects of process parameter variations such as the CO2 amount in the exhaust gas and the feed rate. The presented model is a valuable tool to improve the microbial PHBV production process. Due to the coupling of CO2 online measurements in the exhaust gas to the biomass production, the model has the potential to predict the composition and the current yield of PHBV in the ongoing process.

**Keywords:** bioplastic; copolymerization; polyhydroxyalkanoate; kinetic modeling

#### **1. Introduction**

One suitable alternative to conventional petroleum-based plastics is that of the group of polyhydroxyalkanoates (PHAs) [1,2]. These polyesters stand out because of their favorable processing properties, e.g., their melting behavior or different blending options [3]. Furthermore, these are produced microbially by many bacteria and some archaea using a wide variety of non-fossil carbon sources [1]. There is a repertoire of possible and cheap substrates, such as those of inexpensive sugars in waste streams from the manufacturing industry (juice production, sugar cane processing), volatile fatty acids (VFAs) from biogas plants and wastewater in sewage treatment plants and even CO2 [4–11]. In addition to the diverse possibilities of microbially producing bio-based PHAs, this plastic raw material has another important property: PHAs are degradable under environmental conditions [12]. However, from an economic point of view, the industrial production of PHAs is about five times more expensive than the production of petroleum-based polymers [13]. In addition to an improved extraction and processing of the polymers, a large part of the costs can be saved through optimized bioprocesses with increased PHA yield. This can be achieved by the incorporation of the sophisticated experimental investigation of different process modes or the optimization of substrates and feeding strategies with mathematical modeling [2,14,15]. Furthermore, model approaches represent the basic component in the development of advanced process control and intensification strategies [16,17].

In the research area of PHA production, a large number of models can be found that greatly differ in terms of modeling approach and complexity. Due to the complexity and variability of the bioprocess, there is no universal tool for predicting product yields

**Citation:** Duvigneau, S.; Dürr, R.; Behrens, J.; Kienle, A. Advanced Kinetic Modeling of Bio-co-polymer Poly(3-hydroxybutyrate-*co*-3 hydroxyvalerate) Production Using Fructose and Propionate as Carbon Sources. *Processes* **2021**, *9*, 1260. https://doi.org/10.3390/pr9081260

Academic Editors: Philippe Bogaerts and Alain Vande Wouwer

Received: 29 June 2021 Accepted: 15 July 2021 Published: 21 July 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

regardless of the producing organism, bioreactor or process conditions such as temperature or pH [18]. Some approaches appear promising due to their simplicity and are able to reproduce the concentration curves in a qualitative manner, while they contain only little metabolic information [19–21]. Other approaches take the metabolism into account, however, due to their complexity, they can only be used to a limited extent for model-based process control intensification and are difficult to transfer to other PHA producers [22–27].

Many of the modeling approaches focus on the microbial production of the bestknown representative from the group of PHAs: poly(3-hydroxybutyrate) (PHB). From an industrial point of view, however, the experimentally well-investigated bio-co-polymer poly(3-hydroxybutyrate-*co*-3-hydroxyvalerate) (PHBV) is more interesting, because of its lower melting temperature, higher elongation-to-break values and higher biocompatibility in comparison to PHB [28]. However, only a few model approaches already exist to investigate microbial PHBV production [21,29]. PHBV is also the target product of the present work. In order to develop a universal simulation tool, the mathematical model must contain a balanced amount of metabolic information. Such a modeling approach is rarely found in the literature: in [21], the description of the metabolism was reduced to central points with respect to PHA production (e.g., acetyl-CoA production), which can be found in many organisms in mixed cultures.

In the mathematical model presented here, we applied a time-dependent, kinetic parameter for the formation of residual biomass from fructose and propionic acid in *C. necator*, in order to map the dynamics of the present metabolic activity without detailed metabolic information. As the researchers at the University of Antioquia (Colombia) already have shown [30,31], the online measurement for the CO2 content in exhaust gas serves as an excellent measure of the dynamic growth rate. Here, we also apply this correlation and hence, the model is a suitable candidate for the online estimation of PHBV product yields. In addition, it can be used to predict the polymer composition, since 3 hydroxybutyrate (3HB) and 3-hydroxyvalerate (3HV) amounts in the chains are considered in the model. For the parameter adaption, two data sets from aerobic PHA production in *C. necator* were used, and one experimental setup with only fructose as a carbon source and the other with fructose and propionic acid as carbon sources. In the last part of this manuscript, the model is used in a simulation study to investigate the influence of the feed rate for the propionic acid and of constant CO2 in the exhaust gas on the bio-co-polymer yield and the 3HV proportion in the polymer.

The presented model approach is a suitable candidate for the development of a soft sensor for the online prediction of the polymer yield and composition. This opens new options to increase the flexibility, productivity and quality of the PHBV production process. With further model extensions, e.g., by coupling to polymerization kinetics [32], it will be possible to additionally estimate the chain length distribution during the process to obtain more information about polymer properties.

#### **2. Experimental Methods**

#### *2.1. Mircoorganism and Cultivation Conditions*

*C. necator* (H16, DSM 428) obtained from DSMZ GmbH Braunschweig was used for the fermentations. Bacteria were precultured in a shake flask with 10 vol% LB medium (Carl Roth, Karlsruhe, Germany) at 30 ◦C and 150 rpm. After reaching an optical density of 4 at 600 nm, the bacteria were transferred to an shake flask filled with 10 vol% of M81 medium supplemented with 20 g/L fructose and 1.5 g/L ammonium chloride. The recipe for the Medium 81 can be found in [23] or on the web page of the DSMZ. The M81 preculture was grown until an OD of 4.8 and used as an inoculum for the bioreactors. The fermentations were performed in a DASGIP parallel bioreactor system (Eppendorf AG, Juelich, Germany) with an inoculation OD of 0.4. During the experiments, the pH was kept at 6.8 and the dissolved oxygen (DO) was 70%. The DO measurements were performed with sensors from Mettler Toledo (Gießen, Germany). In the case of fructose as a single carbon source, the pH-control was performed with 2 M H2SO4. During the reactor

experiment with fructose and propionic acid, the pH was stabilized with 20 g/L solution propionic acid as shown in [33]. The detection of the exhaust gas composition was done with the GA4 module of the DASGIP parallel bioreactor system (Eppendorf AG, Juelich, Germany). The initial conditions for the reactor experiments are shown in Table A1. All bioreactor experiments were performed with M81 media at 30 ◦C.

#### *2.2. Determination of Total Biomass*

For the determination of total biomass (TBM), 1 ml culture broth was centrifuged for 10 min at 9600× *g* and 4 ◦C (VWR MicroStar 17R, Pennsylvania, PA, USA). In a second step, the cell pellet was dried over night at 80 ◦C and weighted.

#### *2.3. Enzyme Assay*

By using enzymatic test kits (Kit No. 5390 and No. 10139106035, R-Biopharm AG, Darmstadt, Germany) and following the manufacturer's instructions, ammonium and fructose concentrations were determined from the supernatant of the sample.

#### High-Pressure Liquid Chromatography

Concentrations of 3HB and 3HV were determined by applying the procedure published in [34] using an Agilent 1100 high-performance liquid chromatography (HPLC). For this, 1 mL of the culture broth was alkaline digested as reported in [35]. The samples were filtered through a 0.25 μm nylon membrane and 10 μL were loaded on the reverse phase column (Inertsil 100A ODS-3, 5 μm poresize, 250 × 4.6 mm, MZ-Analysentechnik GmbH, Mainz, Germany) and isocratically eluted with 1 mL·min−<sup>1</sup> at 60 ◦<sup>C</sup> with 92% low concentrated H2SO4 (0.025% solution, Carl Roth, Karlsruhe, Germany) and 8% acetonitrile (Carl Roth, Karlsruhe, Germany). The 3HB and 3HV concentrations in the polymer chains of the samples were determined by using crotonic (Carl Roth, Karlsruhe, Germany) and 2-pentenoic acid standard samples (Sigma Aldrich, St. Louis, MO, USA), respectively. In parallel, a PHBV sample (12% 3HV, Sigma-Aldrich /Merck, Darmstadt, Germany) with known concentration must be measured to calculate the conversion yields *YHB* and *YHB* [34]:

$$\Upsilon\_{HB} = 2 \cdot \frac{c\_{CA}}{c\_{HB}} \, ' \tag{1}$$

$$\Upsilon\_{HV} = 2 \cdot \frac{\mathcal{C}\_{PA}}{\mathcal{C}\_{HV}} \,. \tag{2}$$

Here, the dilution ratio (D) is 2, *cHB* is the known HB and *cHV* the known 3HV concentration of the PHBV test sample. Due to the standard measurement of crotonic acid *cCA* and 2-pentenoic acid *cPA*, the conversion yields *YHB* and *YHV* can be determined, respectively. Detection takes place with a photodiode-array detector (G7115A, Agilent, Waldbronn, Germany) at 210 nm.

#### **3. Kinetic Modeling Approach**

The description of the formation and degradation in the microbial PHA production is an important building block for the complete production process. There are already a number of model candidates for the formation of PHB [18,20,23,36] that can accurately reflect the development of the homopolymer concentration over time. Compared to the homopolymer PHB, the copolymer PHBV has significantly improved processing properties. However, so far, only simple kinetic approaches for the formation of PHBV were developed [21]. Furthermore, the model from Špoljari´c and colleagues was developed for the conversion of fatty acid methyl esters (FAMES) from biofuel to PHBV using lumped metabolic pathways [19]. The model presented here describes the formation and degradation of 3HB and 3HV in the polymer chains using fructose and propionate, two carbon sources that frequently occur in inexpensive residues or can be produced from them, e.g., by using waste streams from juice, cheese and paper production. In our model approach, detailed metabolic reaction pathways were not taken into account to keep the model structure as simple as possible.

The following assumptions were made for the model:


In the following, a set of ordinary differential equations for the dynamics of the system with fructose and propionic acid as substrates and residual biomass, 3HB and 3HV in the polymer chains as products is described. The dynamic state equation for the fructose concentration is given as

$$\begin{aligned} \frac{dc\_{fru}}{dt} &= -k\_1 \cdot b\_{CO\_2}(t) \cdot c\_{res} \cdot c\_{fru} \cdot c\_{\pi} \cdot \text{inh}\_1 \\ &- k\_4 \cdot c\_{res} \cdot c\_{fru} \cdot \text{inh}\_2 \cdot \text{inh}\_3 \\ &- k\_7 \cdot b\_{CO\_2}(t) \cdot c\_{fru} \cdot c\_{res} \\ &- D \cdot c\_{fru} \end{aligned} \tag{3}$$

with:

$$\sinh\_1 = \max\left(0, 1 - \frac{c\_p}{c\_{p,\text{inh}}}\right), \sinh\_2 = \max\left(0, 1 - \frac{c\_n}{c\_n + c\_{n,w}}\right), \sinh\_3 = \max\left(0, 1 - \frac{P\_t}{P\_{l,\text{max}}}\right).$$

Fructose can be metabolized for biomass production with the rate parameter *k*1, the accumulation of 3HB in the polymer with *k*<sup>4</sup> or the conversion to CO2 with *k*7. The growth of biomass through fructose is controlled by the activity coefficient *bCO*<sup>2</sup> (*t*) based on the CO2 ratio in the exhaust gas and inhibited by the concentration of propionate with the term inh1. At a concentration of 1.5 g/L propionic acid (*cp*,*inh*), the substrate uptake for biomass is completely inhibited [33]. Since CO2 in the exhaust gas is often defined as a proportion of the gas composition, we chose the relative CO2 proportion to describe the metabolic activity *bCO*<sup>2</sup> (*t*) as follows:

$$h\_{\rm CO\_2}(t) = \frac{\rm CO\_{2,out}(t)}{\rm CO\_{2,in}} \,. \tag{4}$$

The metabolic activity is described by the quotient of CO2,*out* in the exhaust gas and CO2,*in* in the fresh inlet air. Since *C. necator* is a PHA producer of group 2 according to Novak et al. [18], the build-up of 3HB from fructose begins when there is still a small amount of ammonium in the medium. This effect is modeled by the term inh2. As described in [36], steric effects at high polymer concentrations inhibit the conversion of substrates to PHA (term inh3). According to literature values [1], the maximum achievable amount *Pt*,*max* is 0.89 (89 % of the total biomass).

The inhibitory steric effect is given as the ratio between overall HA concentration and total biomass concentration:

$$P\_t = \frac{(c\_{hb} + c\_{hv})}{(c\_{hb} + c\_{hv} + c\_{rcs})} \,. \tag{5}$$

Finally, the dilution factor in the fed-batch process:

$$D = \frac{F\_{\rm in}}{V} \, , \tag{6}$$

is the ratio of the feed flow rate *Fin* and reactor volume *V*.

For the computational study, a volume balance is necessary:

$$\frac{dV}{dt} = F\_{in}.\tag{7}$$

As for fructose, a state equation can be set up for propionate dynamics:

$$\begin{split} \frac{d\mathbf{c}\_p}{dt} &= -k\_2 \cdot b\_{CO\_2}(t) \cdot \mathbf{c}\_{res} \cdot \mathbf{c}\_p \cdot \mathbf{c}\_n \cdot \mathbf{in} \mathbf{h}\_1 \\ &- (k\_5 + k\_6) \cdot \mathbf{c}\_{res} \cdot \mathbf{c}\_p \cdot \mathbf{in} \mathbf{h}\_2 \cdot \mathbf{in} \mathbf{h}\_3 \\ &- k\_8 \cdot b\_{CO\_2}(t) \cdot \mathbf{c}\_p \cdot \mathbf{c}\_{res} \\ &+ D \cdot \left(\mathbf{c}\_{p,in} - \mathbf{c}\_p\right) .\end{split} \tag{8}$$

It describes the consumption of propionate for biomass with a rate coefficient *k*2, CO2 with *k*<sup>8</sup> and 3HB production with *k*5. In addition to the generation of 3HB, propionate can also be converted into 3HV (*k*8). In fed-batch mode, a propionate solution according to Table A1 is fed to the system with the feed flow rate *Fin*.

For growth, organisms need ammonium. The state equation for the ammonium dynamics is:

$$\begin{split} \frac{dc\_{\rm n}}{dt} &= -c\_{\rm res} \cdot c\_{\rm n} \cdot b\_{\rm CO\_2}(t) \cdot \left(k\_1 \cdot c\_{fru} + k\_2 \cdot c\_p\right) \cdot \text{inh}\_1 \\ &- k\_3 \cdot c\_{\rm res} \cdot c\_{\rm n} \cdot \left(c\_{hb} + c\_{hv}\right) \\ &- D \cdot c\_{\rm n} \cdot \text{.} \end{split} \tag{9}$$

In addition to the ammonium uptake for biomass growth by consuming external carbon sources (first term in Equation (9)), ammonium is needed to convert the biopolymer to catalytically active biomass with the degradation rate parameter *k*3.

The dynamical behavior of residual (non-PHA, catalytically active) biomass is described as follows:

$$\begin{split} \frac{d\boldsymbol{c}\_{\rm res}}{dt} &= \boldsymbol{c}\_{\rm res} \cdot \boldsymbol{c}\_{n} \cdot \left[ \mathbf{in} \mathbf{h}\_{1} \cdot \left( \boldsymbol{k}\_{1} \cdot \boldsymbol{c}\_{f\rm tr} + \boldsymbol{k}\_{2} \cdot \boldsymbol{c}\_{p} \right) \cdot \boldsymbol{b}\_{\rm CO\_{2}} + \boldsymbol{k}\_{3} \cdot \left( \boldsymbol{c}\_{hb} + \boldsymbol{c}\_{h\rm r} \right) \right] \\ &- D \cdot \boldsymbol{c}\_{\rm res} \, . \end{split} \tag{10}$$

Residual biomass is produced through the consumption of external carbon sources such as fructose and propionate and the conversion of 3HB and 3HV from the polymer chains in the presence of ammonium.

The following ordinary differential equations (ODEs) account for the dynamics of the monomers 3HB and 3HV in the polymer chains:

$$\begin{split} \frac{dc\_{hb}}{dt} &= c\_{res} \cdot \text{inh}\_2 \cdot \text{inh}\_3 \cdot \left( k\_4 \cdot c\_{fru} + k\_5 \cdot c\_p \right) \\ &- k\_3 \cdot c\_{res} \cdot c\_u \cdot c\_{hb} - D \cdot c\_{hb} \cdot \text{inh} \end{split} \tag{11}$$

$$\begin{split} \frac{dc\_{hv}}{dt} &= k\_6 \cdot c\_{res} \cdot c\_p \cdot \text{inh}\_2 \cdot \text{inh}\_3 \\ &- k\_3 \cdot c\_{res} \cdot c\_n \cdot c\_{hv} - D \cdot c\_{hv} \,. \end{split} \tag{12}$$

For the accumulation and breakdown of the biopolymer (3HB and 3HV), CO2 formation is negligible, since the metabolic reaction pathways produce only little CO2 compared to the breakdown of sugars and organic acids into catalytically active components of the total biomass.

#### *3.1. Numerical Solution*

3.1.1. Interpolation for Volume, CO2,*out* and Feed Rate

For the simulation of the model, the temporal evolution of the CO2 amount in the exhaust gas from reactor experiments is required. Since the online measurement often fluctuates and does not provide a smooth curve, the data were interpolated to integrate them into the ODE model. For this, a smoothing spline interpolation was carried out with the *MATLAB* command *csaps*. For both experiments, different smoothing factors were evaluated. A smoothing factor of 0.2 was selected for the experiment with fructose as the only carbon source (data set 1), while the data from the reactor experiment with fructose and propionic acid as carbon sources (data set 2) achieved a smooth and well-fitted curve with a smoothing factor of 0.02. For the interpolation, the splines are evaluated at the sampling points. The evaluation was carried out with the *MATLAB* command *ppval*. The curves and online data are shown in Figure 1 for both data sets.

**Figure 1.** Exhaust CO2 for data set 1 (fructose as carbon source, **left** panel, +) and data set 2 (fructose and propionic acid as carbon sources, **right** panel, +) and the interpolations (solid lines).

The feeding of the odd carbon source propionic acid was achieved by pH control. If the pH increases, a certain amount of propionic acid with 20 g/L in the feed is added to the bioreactor to stabilize the pH at 6.8. The pre-implemented PI controller of the DASGIP parallel bioreactor system (Eppendorf AG, Jülich, Germany) was used for this purpose. As for the activity factor, frequent fluctuations are observed because of the special pHdependent feeding strategy and thus the feed rate for propionic acid was interpolated and evaluated in the same way as the CO2 amount in the exhaust gas (Figure 2). Here, a factor of 0.02 delivered a smooth curve.

Furthermore, an interpolation of the volume was necessary for the fed-batch experiment with fructose and propionic acid as carbon sources (data set 2). For this purpose, a polynomial of order 10 was determined with the *MATLAB* command *polyfit* and evaluated with *polyval* at the sampling times. As seen in Figure 3, volume reduction by sampling was also taken into account. Hence, a decrease in volume was recorded despite an average feed rate of approximately 20 mL/h between 15 and 25 h (see Figure 2). In the experiment with fructose as the single carbon source, the reactor was operated in batch mode. Since it was

assumed that the system was ideally mixed, the changes in concentration were only caused by internal sinks and sources (no substrate was pumped in), and the volume reduction due to sampling can be neglected in data set 1. All simulations, interpolations and evaluations were carried out with *MATLAB 2019b*.

**Figure 2.** Experimental feed rate and polynomial for the propionic acid inlet of data set 2 (fructose and propionic acid as carbon sources).

**Figure 3.** Experimental volume (small blue line) and polynomial (black curve) of data set 2 (fructose and propionic acid as carbon sources).

#### 3.1.2. Parameter Identification

For parameter identification, the following objective function was minimized:

$$ESS = \sum\_{i=1}^{n} \left( \frac{\mathbf{x\_{exp}}(t\_i) - \mathbf{x\_{sinh}}(t\_i)}{\max(\mathbf{x\_{exp}})} \right)^2. \tag{13}$$

Here, the error between the simulated **xsim** and experimental data **xexp** at time point *ti* is determined and weighted with the maximum value in the experimental data set.

The kinetic parameters were determined using the algorithm *fmincon* in *MATLAB 2019b*. The ODEs were numerically solved with the algorithm *ode15s* with a relative tolerance of 10−9. To prevent a sub-optimal initial parameter set a multi-start approach with N = 10,000 was applied. To further validate the resulting parameter set obtained by the local optimization strategy, the global optimization algorithm differential evolution (DE) was selected [38]. The resulting set of parameters can be found in Table A1.

To prove the parameter identifiability, profile likelihoods are determined for the parameter set [39]. Appendix B shows all profile likelihoods and the corresponding likelihood-based confidence intervals. All profile likelihoods show a distinct minimum at the estimated parameter set and thus all parameters and the model itself are (locally) identifiable. The confidence intervals are obtained by calculating the value of a *χ*2-distribution with a confidence level *α* = 0.95 and one degree of freedom as proposed in [39].

#### **4. Results**

#### *4.1. Identification Using Different Data Sets*

The kinetic model was adapted to two data sets from bioreactor experiments with a working volume of 1.2 L. In the first data set, fructose was the only carbon source that was metabolized under aerobic conditions. The second data set was also obtained under aerobic conditions with fructose and propionic acid as carbon sources. Here, propionic acid was added via a pH-regulated feed as proposed by Kim and coworkers [33] aiming for a constant pH value of 6.8. From the available data of these sets, the online data for the CO2 content in the inflow and in the exhaust gas, the feed rate for the propionic acid and the exact volume considering the sampling volume were used in the case of data set 2. In the case of the data set 1, online data for the CO2 content in the inflow and in the exhaust gas were also used, but a constant volume and a batch mode (F*in* =0) were considered. The CO2 content in the exhaust gas, the feed stream for the propionic acid and the reactor volume were approximated as described in Section 3.1.1. The smoothed measurement data were used for the model simulation. Furthermore, the concentrations for total biomass, biopolymer, fructose and propionic acid were determined offline in both data sets (see Section 2).

The dynamic behavior for the conversion of the substrates from data set 1 (only fructose) can be reproduced well with the present model (solid lines, Figure 4a). The model shows larger deviations for the substrates from data set 2 (fructose and propionic acid), especially in the last time segment from 25 h (dashed lines, Figure 4b). On the one hand, this is due to the approximation of the inflow rate for propionic acid (see Figure 2), and on the other hand, the measurement of the propionic acid in the medium becomes more difficult. It seems to be, that there are more and more apoptosis fragments, e.g., matrix, RNA and proteins in the culture supernatant that disrupt the signal obtained by HPLC (own experimental findings). These fragments could be avoided by elaborate sample preparation before HPLC measurement, e.g., additional filtration, boiling procedures or the supplementation of organic acids. Furthermore, the chromatographic peak consisting of propionic acid can be be fractionated and separated from impurities by a second HPLC run with an adjusted mobile phase.

The model for the case of fructose as a single substrate can reproduce the production and depletion of total biomass and 3HB with sufficient accuracy (Figure 5a). The same applies to the case with fructose and propionic acid as carbon sources (Figure 5b). In particular, the degradation of 3HB and 3HV in the polymer chains after an NH4Cl shot at 24 h can be mapped very well by the model. This property is important when working with waste streams which also contain nitrogen sources and which are added during the ongoing process in order to keep the carbon in excess. Since the model can reproduce the product concentrations very well, it can further be used for a simulation study.

(**a**) *Data set 1 (fructose only)—substrates* (**b**) *Data set 2 (fructose + propionic acid)—substrates*

**Figure 4.** Consumption of substrates for two feeding scenarios: (**a**) fructose as a single carbon source without additional feeding (batch); (**b**) fructose and propionic acid as carbon sources with pH-dependent propionic acid feeding.

(**a**) *Data set 1 (fructose only) - products* (**b**) *Data set 2 (fructose + propionic acid) - Products*

**Figure 5.** Production and degradation of polyhydroxyalkanoate (PHA) monomers (chb, chv) and total biomass (cbio = cres + chb + chv).

#### *4.2. Computational Study*

In the following, the identified model was used to investigate the maximum product concentrations by applying different constant CO2 and propionic acid feed profiles. A constant feed rate can easily be implemented on every bioreactor system with continuous pumps. The constant CO2 in the exhaust gas represents a complex control task, as CO2 is produced by the bacteria themselves (autogenous CO2). However, the research work in [40] shows the feasibility of this control task which can also be used in future validation experiments.

In the simulation study (Figure 6), a constant CO2 proportion in the exhaust gas and a constant feed rate for propionate were assumed for the process time. The propionic acid concentration in the feed was set to 20 g/L as in data set 2. Figure 6a shows the maximum total 3HB and 3HV in the polymer chains P*t*\* as a function of CO2 in the exhaust gas and the feed rate for propionic acid. The inhibitory area is clearly visible on the left side. In this range, the propionic acid concentration in the medium becomes too high so that growth is inhibited. In the case of increased CO2 values in the exhaust gas, the feed rate can also be increased without triggering growth inhibition. This behavior can be justified as follows: a higher exhaust gas value for CO2 produced by the microorganisms (autogenous CO2) is triggered by an increased uptake of substrates from the medium. In *C. necator*, the degradation of pentoses and hexoses takes place via the Entner–Doudoroff (ED) pathway. Furthermore, the tricarboxylic acid cycle (TCA) is mainly responsible for the generation of energy and precursor molecules for biomass synthesis from the precursors of metabolic sugar degradation and organic acids. Both the ED pathway and the TCA generate CO2 as a by-product, which can be found in the exhaust gas values of the simulation study. With other words, increased autogenous CO2 in the exhaust gas can be translated into stronger residual biomass growth with increased substrate uptake. As a result, the growth inhibition with increased CO2 in the exhaust gas only occurs at higher feeding rates for propionic acid. In addition to the beneficial effect of higher exhaust CO2, the total biopolymer concentration decreases in the non-inhibitory area. This effect occurs because the CO2 output is related to the higher residual biomass growth and hence, the substrates were less translated into biopolymers.

(**a**) *Best total 3-hydroxybutyrate (3HB) and 3-hydroxyvalerate (3HV) concentration*

(**b**) *Maximum possible 3HV concentration after 60 h*

**Figure 6.** Simulation study with constant values for CO2 and F*in*.

The highest total biopolymer concentration of 12.5 g/L within 60 h simulation time can be achieved without feeding propionic acid to the system and with a very low CO2 content in the exhaust gas. However, without an active feed rate, no 3HV will be produced and the 3HV content is decisive for improved processing of copolymers compared to their homopolymers. The maximum HV concentration in the polymer is shown in Figure 6b for different constant autogenous CO2 and propionic acid feeding rates. Here, the inhibiting region can be seen again due to a high propionate concentration in the medium. For a high 3HV concentration, our model predicts feed rates between 12 and 40 mL/h depending on the CO2 in the exhaust. In general, less than 2.5 % CO2 in the exhaust gas leads to higher 3HV concentrations. As our simulation results show, the feed rate for propionic acid must be chosen very carefully because of their strong correlation with the CO2 value: overly high feed rates at higher CO2 in the exhaust gas lead to less residual biomass and a decreased 3HV concentration (see Figure A1). Further characteristic values of the simulation study are shown in Figure A1 (e.g., residual biomass at the maximum total 3HB and 3HV concentration, total monomer/total biomass ratio).

Three exemplary time courses for different production goals were illustrated in Figure 7. In case A, the dynamic behavior to achieve a high total polymer concentration (3HB + 3HV) with the given initial conditions is shown. For this purpose, the feed rate was set to 0 mL/h and the CO2 amount to 1 %. For case B, the same CO2 value as in case A was applied together with a feed rate of 25 mL/h to show an example time course for a high 3HV concentration. In addition to the two preferred fermentation results (cases A and B, Figure 7), the inhibitory case was also shown (case C, Figure 7) by increasing the feed rate to 105 mL/h.

**Figure 7.** Three example cases of dynamical behavior applying constant exhaust CO2 and feed rate for propionic acid (profile study 1). Legend: A, maximum total biopolymer concentration (chb + chv); B, high 3HV concentration (chv); and C, inhibition caused by propionic acid (cp).

#### **5. Concluding Remarks**

In this manuscript, a model approach is presented which enables the integration of online data for the estimation of the yield and the composition of the copolymer PHBV in *C. necator*. Compared to other approaches, no genome-scale metabolic networks or reduced variants are necessary [23,41], since it is a pure kinetic approach that describes changes in the metabolism due to a CO2-dependent biomass production rate that changes over time. Despite the lack of detailed metabolic information, our kinetic model can display the data sets with fructose as a single substrate and fructose and propionate as substrate with high accuracy. The model approach has similar complexity as the models presented in Koller et al. [36] and Dias et al. [21], whereby the metabolic activity that controls the biomass growth can be described by the CO2 profile in the exhaust gas. Because of the compact model structure and smart coupling to available online measurements the approach can be applied as a soft sensor for the prediction of biomass, substrates, product yield and the composition of the biopolymer (e.g., 3HV/3HB ratio). First steps in this direction where already successfully applied for a vinasse-molasse PHA process [30,31]. In comparison to the model approach used in [30,31], our more complex model represents an excellent basis to predict 3HV and 3HB amounts in the polymer chains. In addition to the application options during the process, CO2 profiles can be estimated in optimization studies with predefined concentrations profiles for the substrates and products. As a first step, we investigated in our simulation study the effect of constant CO2 fractions in the exhaust gas and constant feed rates for propionic acid on biopolymer yield and composition. The lower the feed rate was set and the less CO2 was in the exhaust gas, the more PHA was produced but with less 3HV content in the polymer. A suitable feed rate for propionic acid input was predicted to be between 12 and 40 mL/h in order to achieve high 3HV concentration in the final co-biopolymer. Furthermore, the model approach can be used for the design of observers and state estimators for the reconstruction of non-measurable states. A first example in this direction was recently presented by Carius and coworkers [16]. Here, an unscented Kalman filter and a moving horizon estimation based on a hybrid cybernetic PHB model [23] was designed and evaluated.

Future work should focus on the design of model-based soft sensor approaches, since this will enable the reliable online estimation of the PHBV content using measurements of the exhaust gas online without the need for additional expensive hardware sensors. Furthermore, the transferability of the kinetic model to other PHA producers must be researched. Here, the focus should be on PHA producers, which are already producing PHA under growth conditions, e.g., *C. necator* DSM 515, as the model was designed for this group of bacteria. Furthermore, it should be checked whether a characteristic CO2 profile occurs during the accumulation of other copolymer building blocks, such as 4 hydroxybutyrate. Finally, a control concept should be developed which is able to keep CO2 in the exhaust gas on a desired level over a longer period of time. The work of Shang and colleagues [40] shows that it is in principle possible to control a process parameter that is strongly influenced or caused by the bacteria. The CO2 amount in the exhaust gas was adjusted in the work of Shang and coworkers in order to investigate the inhibitory effects of CO2. Such an effect has not yet been taken into account in the model presented in this manuscript, but should be considered in future model extensions.

Overall, our model approach provides the basis for a broad range of possible future applications and will help make the production process of biopolymers more reliable and less expensive.

**Author Contributions:** S.D. and R.D. conceived and designed research. S.D. and J.B. conducted the experiments. S.D. analyzed and visualized the data and wrote the original draft of the manuscript. R.D. and A.K. supervised the work and edited the original draft. All authors have read and agreed to the published version of the manuscript.

**Funding:** We would like to acknowledge the EU-programme ERDF (European Regional Development Fund) for the funding of the project DIGIPOL.

**Institutional Review Board Statement:** Not applicable

**Informed Consent Statement:** Not applicable

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

**Appendix A. Additional Figures—Profile Study**

(**a**)*3HV/total biopolymer ratio at maximum total biopolymer concentration after 60 h*

(**b**)*Residual biomass concentration at highest total biopolymer concentration after 60 h*

(**c**)*Biopolymer concentration/total biomass ratio at maximum total biopolymer concentration after 60 h*

(**d**)*Time where the maximum biopolymer value occurs*

**Figure A1.** Simulation study with constant values for CO2 and F*in*—3HV content (**a**), residual biomass (**b**), PHA/BTM ratio (**c**), time (**d**).

**Figure A2.** Reactor volume V, feed rate F*in*, dilution rate D and CO2 in the exhaust for the three example cases. Legend: A, maximum total biopolymer concentration; B, high 3HV concentration; and C, inhibition caused by propionic acid.

**Appendix B. Parameter Identifiability**

**Figure A3.** Profile likelihood of k1. Horizontal dashed line represents the explained sum of squares (ESS) value of the 95% confidence interval.

**Figure A4.** Profile likelihood of k2. Horizontal dashed line represents the ESS value of the 95% confidence interval.

**Figure A5.** Profile likelihood of k3. Horizontal dashed line represents the ESS value of the 95% confidence interval.

**Figure A6.** Profile likelihood of k4. Horizontal dashed line represents the ESS value of the 95% confidence interval.

**Figure A7.** Profile likelihood of k5. Horizontal dashed line represents the ESS value of the 95% confidence interval.

**Figure A8.** Profile likelihood of k6. Horizontal dashed line represents the ESS value of the 95% confidence interval.

**Figure A9.** Profile likelihood of k7. Horizontal dashed line represents the ESS value of the 95% confidence interval.

**Figure A10.** Profile likelihood of k8. Horizontal dashed line represents the ESS value of the 95% confidence interval.

#### **Appendix C. Parameter Values and Initial Conditions**

**Table A1.** Kinetic parameters, states and variables.



**Table A1.** *Cont.*

#### **References**

