Multiscale dynamic modeling and simulation of a biorefinery

Abstract A biorefinery comprises a variety of process steps to synthesize products from sustainable natural resources. Dynamic plant‐wide simulation enhances the process understanding, leads to improved cost efficiency and enables model‐based operation and control. It is thereby important for an increased competitiveness to conventional processes. To this end, we developed a Modelica library with replaceable building blocks that allow dynamic modeling of an entire biorefinery. For the microbial conversion step, we built on the dynamic flux balance analysis (DFBA) approach to formulate process models for the simulation of cellular metabolism under changing environmental conditions. The resulting system of differential‐algebraic equations with embedded optimization criteria (DAEO) is solved by a tailor‐made toolbox. In summary, our modeling framework comprises three major pillars: A Modelica library of dynamic unit operations, an easy‐to‐use interface to formulate DFBA process models and a DAEO toolbox that allows simulation with standard environments based on the Modelica modeling language. A biorefinery model for dynamic simulation of the OrganoCat pretreatment process and microbial conversion of the resulting feedstock by Corynebacterium glutamicum serves as case study to demonstrate its practical relevance.

modeling as an important tool at different stages of process development.
During early-stage process design, a large number of product and process alternatives is analyzed and compared with respect to economic and environmental criteria using mathematical optimization. Simple models are used, for example, yield coefficients for reactions (Bao et al., 2011;Voll & Marquardt, 2012) and short-cut methods for separation steps (Ulonska, Skiborowski, Mitsos, & Viell, 2016), combined with detailed models (e.g., equilibrium and kinetic models) where data are available (Kelloway & Daoutidis, 2014;Martín & Grossmann, 2013) to identify promising feedstocks and reaction pathways. More detailed, steady-state, plant-wide process models are used to investigate and improve individual biorefinery concepts. For example, Tay, Kheireddine, Ng, El-Halwagi, and Tan (2011) investigated gasification of lignocellulosic biomass in an integrated biorefinery using a superstructure optimization framework with models of different complexity.
In the next step, these biorefinery concepts are experimentally investigated in pilot plants requiring dynamic process models to improve process understanding and to allow model-based process operation and control. As the number of considered process alternatives decreases, more detailed models can be formulated.
While this is true for most process steps, the complexity of microbial transformations is usually not represented on the same level of detail. In the past, mostly unstructured (so called black-box) models were used to approximate the reaction kinetics of microbial transformations.
However, this class of models is not suitable to describe and predict the complex intracellular metabolism of microorganisms cultivated in bioreactors for the purpose of bio-based production. Clearly, to understand these multilevel interactions in a real quantitative manner, mechanistic pathway modeling in combination with multi-omics analytics would be required (Wiechert & Noack, 2011). Although a couple of promising examples in this direction do exist (Gonçalves et al., 2017;Hameri, Fengos, Ataman, Miskovic, & Hatzimanikatis, 2018;Noack, Voges, Gätgens, & Wiechert, 2017;Zieringer & Takors, 2018) this approach is still hampered by the availability of in vivo enzyme kinetic data for relevant metabolic reactions.
As an intermediate solution, dynamic flux balance analysis (DFBA) enables the simulation of consecutive, steady-state metabolic flux distributions under changing environmental conditions (Mahadevan, Edwards, & Doyle, 2002;Palsson, 2008;Zhao, Noack, Wiechert, & von Lieres, 2017). DFBA is based on a linear equation system covering the stoichiometry of all intracellular reactions within defined network boundaries and a differential equation system covering the timedependent behavior of extracellular biomass, substrate, and product concentration (Mahadevan et al., 2002). In most cases, the system of equations is under-determined. Formulation of a cellular optimization criterion, for example, maximization of growth or ATP production, allows choosing among the set of solutions to determine one particular metabolic flux distribution. Finally, the DFBA model represents a system of ordinary differential equations with embedded optimization criteria (ODEO; Zhao et al., 2017).
Unfortunately, this type of mathematical problem is hard to solve and requires specialized solution methods that are commonly categorized as static optimization approach (SOA), dynamic optimization approach (DOA), and direct approach (DA). SOA divides the time horizon into intervals, solves the embedded optimization problem at the beginning of each interval and performs integration assuming a fixed flux distribution for the respective interval (Mahadevan et al., 2002). DOA discretizes the dynamic equations using collocation methods to obtain a nonlinear problem (NLP) that needs to be solved only once (Mahadevan et al., 2002). DA uses a solver to evaluate the embedded optimization problem in each time step during integration (Hanly & Henson, 2011;Hjersted & Henson, 2006;Zhuang et al., 2011). Höffner, Harwood, andBarton (2013) have implemented a DFBA simulator that uses first-order optimality conditions to reformulate the embedded optimization problem into a system of differential-algebraic equations (DAEs).
Based on these solution techniques, the DFBA approach was successfully used for various applications, for example, for the investigation of the diauxic growth of Escherichia coli on D-glucose (Mahadevan et al., 2002), the optimization of a fed-batch fermentation process with Saccharomyces cerevisiae on D-glucose (Hjersted & Henson, 2006), and the optimization of co-cultures with both model organisms on mixtures of D-glucose and D-xylose (Hanly & Henson, 2011).
One particular advantage of the DFBA approach is its ability to account for changing environmental conditions during (fed-)batch operation. Therefore, we chose the multiscale modeling approach by embedding a DFBA model into a plantwide biorefinery simulation.
Mathematically, this leads to a system of differential-algebraic equations with embedded optimization criteria (DAEO).
In this study, we present a framework that enables dynamic We implemented a library of dynamic models for unit operations that are widely used in chemical processes, for example, flash unit, decanter, and distillation column, using the object-oriented modeling language Modelica (Fritzson, 2015). Modelica is well-suited for modeling biorefineries as it allows replacement of simple models with more detailed ones as new insights become available and also enables replacement of entire process steps if desired. The library uses the components Modelica.Fluid and Modelica.Media of the Modelica Standard Library (Casella, Otter, Proelss, Richter, & Tummescheit, 2006;Elmqvist, Tummescheit, & Otter, 2003). The former provides an interface for one-and zero-dimensional physical modeling of thermohydraulic components (Casella et al., 2006). The component equations (e.g., mass and energy balance equations) are decoupled from fluid property equations (e.g., calculation of specific enthalpy or density). In this way, the same component model can be used with different fluids (Casella et al., 2006). The Modelica.Media library defines the fluid and allows implementation of different models such as ideal gases and multiphase systems (Elmqvist et al., 2003).
In  (Biegler, 2010;Gopal & Biegler, 1999;Sahlodin, Watson, & Barton, 2016) and LLE (Müller & Marquardt, 1997;Ploch, Glass, Bremen, Hannemann-Tamás, & Mitsos, 2019). The phase calculation may again be performed externally using a suitable package. Reactors were modeled with input-output relations based on experimental data or kinetic models if available. The piping between two building blocks was not considered, that is, energy and pressure loss were neglected. The graphical interface allows generation of process models at the flowsheet level.

| DFBA process model for microbial transformation
The dynamic model for the microbial conversion step is based on the DFBA approach given by where x are the n x differential states describing the extracellular environment with initial conditions given by x 0 . The differential equations f depend on the optimal solution v of an embedded optimization problem. Note that a typical DFBA process model does not include any algebraic equations outside the embedded optimization problem. Therefore, we refer to equation system (1a) as ODEO.
The optimization problem is based on classical FBA. The stoichiometry of the reaction network is defined by the stoichiometric matrix  (Hanly & Henson, 2011;Varma & Palsson, 1994) or ATP production (Raghunathan, Pérez-Correa, Agosin, & Biegler, 2006) or nonlinear, for example, minimization of the overall intracellular flux (Schuetz, Kuepfer, & Sauer, 2007) or maximization of biomass yield per flux unit (Zhao et al., 2017). In this study, we used a linear objective function that can be defined by (ˆ) = −ĥ v cv T , where c is the cost vector.

| Framework for plant-wide biorefinery modeling and simulation
We established a framework for dynamic modeling and simulation of biorefineries that covers several steps, which are illustrated in Figure 1 and briefly described in the following: 1. We implemented a model library of building blocks in Modelica enabling dynamic modeling of early process steps such as pretreatment of lignocellulosic biomass (see Methods section for more details).
F I G U R E 1 Framework for plant-wide biorefinery modeling and simulation. For modeling biomass pretreatment processes a library of dynamic unit operations was implemented in Modelica. Thermodynamic property parameters were extracted from Aspen Plus and converted to a Modelica record. For modeling the microbial transformation step the visualization tool Omix is used and the resulting metabolic network is parsed to Modelica code. An external program interprets the code and creates the DFBA submodel. The external program applies the direct solution method via an interface to LP solver CPLEX. This way, a simple and easy-to-use interface is obtained and Dymola or OpenModelica can be used to solve the resulting multiscale biorefinery model (DAEO In the following we applied this framework to simulate the dynamics of the OrganoCat pretreatment process and the subsequent conversion of the resulting sugar streams by Corynebacterium glutamicum in different scenarios.

| Dynamic model of OrganoCat process
Plant-derived biomass as feedstock has to be pretreated to make the native biomolecules accessible for subsequent chemical or microbial transformation. The effectiveness of a pretreatment depends on many factors such as type of biomass, process conditions, formation of unwanted degradation products, recyclability of catalysts, energy and catalyst costs, and many more. An overview of various pretreatment processes for lignocellulosic biomass, their advantages and drawbacks as well as an economic assessment can be found elsewhere (Hendriks & Zeeman, 2009;Mosier et al., 2005).
In this study, we considered the OrganoCat process (Grande et al., 2015;Viell, Harwardt, Seiler, & Marquardt, 2013;Vom Stein et al., 2011) as an example process for pretreatment of lignocellulosic biomass. It uses a biphasic solvent system for selective depolymerization of lignocellulosic biomass into three separate process streams under mild reaction conditions (Vom Stein et al., 2011). Two of these process streams consist of a mixture of C5 (mainly D-xylose) and C6 (mainly D-glucose) sugars with defined ratios that are envisioned to be processed in microbial transformation steps to yield high-value products. To discuss the impact of different sugar ratios we define the fraction of total carbon atoms stemming from D-xylose as with c C5 and c C6 denoting the concentrations of D-xylose and D-glucose, respectively. Note that = 0 C5 ϕ refers to a pure D-glucose The dynamic model of the repetitive-batch OrganoCat process that was formulated in this study is based on the conceptual process design of Viell et al. (2013) and the process variation proposed by Grande et al. (2015). Thermodynamic parameters were equal to the ones used in steady-state simulation (Viell et al., 2013) except for the LLE between water and 2-MTHF that was described by the NRTL model using up-to-date parameter values from literature (Glass, Aigner, Viell, Jupke, & Mitsos, 2017). By applying ideal dynamic unit operations and some additional assumptions (instantaneous flash units, residence time reactors for enzymatic hydrolysis and oxalic acid crystallization) a start-up process was implemented ( Figure 2).
In the first process step, lignocellulosic biomass (e.g., beech wood) is  (Viell et al., 2013). Both the remaining aqueous and organic phase may be re-used in the biphasic reactor.

| Dynamic model of microbial transformation with C. glutamicum
For the microbial transformation step, we considered the platform organism C. glutamicum that is known for its capability to produce a variety of value-added products including amino acids, organic acids, aromatic compounds, and proteins (Baritugo et al., 2018). Wild-type C. glutamicum cannot naturally utilize D-xylose, but in recent years the isomerase pathway and the Weimberg pathway as oxidative strategies for D-xylose metabolization could be functionally implemented into this organism (Meiswinkel, Gopinath, Lindner, Nampoothiri, & Wendisch, 2013;Radek et al., 2014). Interestingly, there are two further alternative routes for D-xylose assimilation known, namely, the oxido-reductase and the DAHMS pathway (see, e.g., Valdehuesa et al., 2018, for a general overview on these pathways), but none of these has been tested for C. glutamicum so far.
F I G U R E 2 Dymola flowsheet of implemented OrganoCat pretreatment process. The biomass is exposed to a biphasic reaction medium (R1) to dissolve lignin and hemicellulosic sugars in the reaction medium. The cellulose fraction remains solid and is subsequently washed (W) and converted into a D-glucose rich sugar stream (GLC) via enzymatic hydrolysis (EH). The liquid outlet of R1 is decanted (DEC) to obtain an organic fraction with dissolved lignin and a D-xylose rich sugar stream (XYL) from the aqueous fraction. Crystallization (CRY) is used to recycle catalyst oxalic acid from the aqueous fraction [Color figure can be viewed at wileyonlinelibrary.com] Therefore we were interested to study each of these pathways The balanced extracellular species are the carbon sources D-glucose (GLC) and D-xylose (XYL), the biomass (X) as well as the model target product succinate (SUCC) and by-products lactate (LAC) and acetate (ACE). Uptake rates for D-glucose and D-xylose were modeled by classical Michaelis-Menten kinetics. In addition, an upper limit for the overall uptake rate of carbon sources was imposed to account for potential transport limitations: Maximization of biomass growth via the specific growth rate μ was used as optimization criterion in all simulations. The full DFBA model including kinetic parameter values is given in Appendix A.

| Case study: Growth of C. glutamicum on OrganoCat media
In a first simulation study, we examined the potential of our , where t B is the final batch time defined by depletion of F I G U R E 3 Simulation results for dynamic pretreatment process (repetitive-batch OrganoCat). Hold-ups of the respective unit operations (as indicated in Figure 2) are shown. The reactor R1 reuses the same medium three times before it is directed to a decanter (DEC). In between each cycle, the solid cellulose pulp is removed. It is further converted to monomeric sugars via enzymatic hydrolysis (EH) yielding a sugar stream with high D-glucose content (GLC, = 0.09  In our third simulation study, we were interested in an optimal production scenario for succinate under microaerobic fed-batch conditions. Therefore we coupled the optimal media design for biomass growth in a first batch phase (i.e., = 0.2 C5 ϕ , cf. Figure  The highest product titer, however, comes at low selectivity (amount of succinate per total amount of carbonic acids formed) due to simultaneous formation of acetate. Consequently, based on the assumed strain design (i.e., C. glutamicum wild type harboring four potentially active D-xylose assimilation pathways, but no further strain modification, Figure 4) microaerobic production of succinate always leads to a trade-off between high product titer and high selectivity.
The combination of optimal batch and fed-batch operation with highest product titer is shown on the right-hand side of Figure 7.
During the batch phase ( ⩽ ⩽ t 0 11 hr), a high growth rate is obtained by metabolization of D-glucose at the highest possible rate. In

| CONCLUSION AND OUTLOOK
We presented a multiscale and multidisciplinary modeling framework that allows dynamic modeling and simulation of complex biorefinery processes. The developed Modelica library of (bio)chemical unit operations permits easy adaptation to describe different subprocesses for the conversion of renewable feedstock into value-added products or to replace (sub)models when new insights and more detailed models become available. The DFBA approach enables realistic modeling and simulation of microbial transformation steps and accounts for changing environmental conditions during dynamic operation. Our solution strategy for the embedded optimization problem is easy to use and allows problem formulation and modification in Modelica language.
Finally, our case study showed that sugar mixtures from OrganoCat pretreatment are well-suited carbon sources for bio-based production with the platform organism C. glutamicum.
The direct approach is suitable for solving the embedded optimization problem during dynamic simulation of DFBA models.
However, it is important to note that numerical difficulties may arise because active set changes of the embedded optimization problem are hidden from the integrator. In addition, the direct solution approach requires the simulation software to use a numerical approximation technique (e.g., finite differences) to obtain derivative sugars (cf. Equation (2)). Vertical dashed lines indicate the two potential mixtures derived from OrganoCat pretreatment process (cf. Figure 2)  on succinate production with Corynebacterium glutamicum under microaerobic fed-batch conditions in terms of final titer, space-time yield and selectivity (left-hand side) and model prediction for combination of optimal batch and fed-batch operation for highest product titer (right-hand side). The set-up for the simulation study was as follows: Constant volumetric feeding rate of =

| 2571
substitution of some (sub)models of the pretreatment process by more detailed models is necessary for in-depth analysis of the influence between the various process steps and identification of reasonable feedback loops. For instance, a more detailed model of the enzymatic hydrolysis step would enable a feedback loop regarding the required sugar ratios for microbial transformation and appears to be a promising handle to improve the overall process performance. In addition, the presented tool may be used to support scale-up of biorefinery processes, for instance, in terms of plant-wide model-predictive control, sizing of equipment, and identification of potential bottlenecks (e.g., availability of reactants, important recycle streams, accumulation of potential inhibitors for enzymatic hydrolysis, or microbial transformation).

ACKNOWLEDGMENTS
This study was financially supported by the Bioeconomy Science

CONFLICT OF INTERESTS
All authors declare no competing interests.

AUTHOR CONTRIBUTIONS
T. P. implemented the Modelica models and the DAEO toolbox, conducted the case studies, and prepared the manuscript. X. Z. All authors read and approved the submitted manuscript.

DATA AVAILABILITY
The Modelica library, the Omix converter, and the DAEO toolbox will be made available for interested users.