 Research
 Open Access
 Published:
Advanced integration of fluid dynamics and photosynthetic reaction kinetics for microalgae culture systems
BMC Systems Biology volume 12, Article number: 93 (2018)
Abstract
Background
Photosynthetic microalgae have been in the spotlight of biotechnological production (biofuels, lipids, etc), however, current barriers in mass cultivation of microalgae are limiting its successful industrialization. Therefore, a mathematical model integrating both the biological and hydrodynamical parts of the cultivation process may improve our understanding of relevant phenomena, leading to further optimization of the microalgae cultivation.
Results
We introduce a unified multidisciplinary simulation tool for microalgae culture systems, particularly the photobioreactors. Our approach describes changes of cell growth determined by dynamics of heterogeneous environmental conditions such as irradiation and mixing of the culture. Presented framework consists of (i) a simplified model of microalgae growth in a culture system (the advectiondiffusionreaction system within a phenomenological model of photosynthesis and photoinhibition), (ii) the fluid dynamics (NavierStokes equations), and (iii) the irradiance field description (BeerLambert law). To validate the method, a simple case study leading to hydrodynamically induced fluctuating light conditions was chosen. The integration of computational fluid dynamics (ANSYS Fluent) revealed the inner property of the system, the flashing light enhancement phenomenon, known from experiments.
Conclusion
Our physically accurate model of microalgae culture naturally exhibits features of real system, can be applied to any geometry of microalgae mass cultivation and thus is suitable for biotechnological applications.
Background
The global warming, gradual oil fuel depletion, higher demands for energy and food consumptions are some of the current challenges our society is facing. Thus, looking for alternative resources is a hot topic. After the failure of the first generation biofuels based on corn and soya, which created a food shortage in the third world, the scientific community focused on simple photosynthetic organisms, cyanobacteria and microalgae. These organisms have high photosynthetic efficiency, biomass growth and lipids content [1].
Cyanobacteria and microalgae are very versatile microorganisms which have attracted a significant amount of scientific interest in the last decades [2–4]. These unicellular oxygenic organisms absorb photons of specific wavelengths to fix inorganic carbon (CO_{2}). Carbon molecules are then utilized in the CalvinBenson cycle, producing complex organic molecules in dark photosynthetic reactions and thus reducing the concentration of carbon dioxide in the environment. In addition, cyanobacteria and microalgae consume inorganic phosphorus and nitrogen which could be used in the processes of wastewater treatment [5]. The current scientific focus is to employ microalgae and cyanobacteria in biotechnology field as a potential biofuel source, for example a genetic modification increasing lipid content [6]. The other applications, e.g., pharmaceutical (bioactive metabolites), nutritional (animal feed and human food and supplements), are less important compared to the renewable energy business.
However, there are technological and knowledgebased barriers preventing optimized mass cultivation of photosynthetic microorganisms. In particular, reliable approaches allowing qualitative predictions for simulation of microbial biomass production in microalgae culture systems (MCS), including the photobioreactors (PBR), are rather slowly emerging than being well established [7–9]. The cause of apparent lack of accurate models of microalgae growth is in the complexity of threedimensional multiphase (liquidgassolid) flow dynamics, irradiance distribution within MCS and multilevel processes of cellular functions. A reliable model has to integrate these phenomena which interact across various timescales. Therefore, it is necessary to solve theoretical (multitimescale phenomena) and practical (high computational requirements) problems. As a result, the complex interactions within the MCS are often integrated in inadequate way, leading to oversimplified models [10]. Especially, the correct integration of a computational fluid dynamics (CFD) code and photosynthetic reaction kinetics is essential for meaningful solution [11]. Finally, current models, describing general MCS or specific PBR performance analysis, have usually no connection to the bioproduction itself, e.g., [4, 12], and references within.
In this study, we present a unified modeling framework of fluid dynamics and photosynthetic reaction kinetics for microalgae culture system of any geometry. The framework is validated by a numerical simulation of microalgae growth in a 2D square cavity with one moving wall, a hypothetical system triggering a vortex flow. This system is able to reproduce the hydrodynamically induced high frequency lightdark cycles regime causing the socalled flashing light enhancement [13–15].
Methods
We present a simulation tool for MCS (photobioreactor, open or raceway ponds) to connect the key aspect of real system, i.e., microalgae cell growth under changing hydrodynamical and optical conditions. This framework consists of (i) a state system – mass balance equations in form of advectiondiffusionreaction partial differential equations – PDEs (13), (ii) the fluid flow equations, i.e., NavierStokes equations (4), and (iii) the irradiance distribution inside MCS according to BeerLambert’s law (5), see [8, 16].
All mentioned parts of the introduced MCS model are interconnected, i.e., the mass balance equations for the state variables (species characteristics, nutrients, gases, etc) are solved simultaneously with the fluid dynamics (momentum balances, continuity equations). Nevertheless, we assume that the stationary flow field inside the MCS is not affected by mass transfer and reactions [17]. This assumption allows the detachment of biological and environmental states and thus biological and environmental parts (models) can be analyzed with different numerical methods and with different spatialtemporal discretizations, dramatically reducing the computational demands. This detachment could be complete (for the stationary flow field in a continuous system) or stepwise, e.g., reflecting some sequence of quasisteady states in a production system operated in batch mode.
State model
The biological part of modeling framework is based on the mass balance equations for the state variables describing the transport and reaction among the species or compounds [18]:
where c_{i}=c_{i}(x,t), is a conservative quantity (concentration or cell density), v is a velocity of flow field ruled by the fluiddynamic model, cf. (4), and x∈Ω⊂R^{3} stands for a position vector in a coordinate system (e.g., Cartesian or cylindrical). The dispersion coefficient D_{e} is a secondorder tensor, which corresponds to the diffusion coefficient in microstructure description. D_{e} is an empirical parameter describing mixing in the system and it is influenced by the molecular diffusion and velocity profile, i.e., it is not strictly a material constant.
The reaction kinetics and temporal changes of the reacting species are described by the reaction term R(c_{i}) and source terms S(c_{i}), respectively. The source term, e.g., the external load of nutrients into MCS, is usually modeled as a corresponding boundary condition. However, in order to (i) simplify the analytic study of the optimal solution existence, see [19], and (ii) respect that the location of discharge of some material could be inside the domain Ω, we employed the particular form of (1). The initial condition and boundary condition (impermeability of the domain boundary, e.g., PBR walls) are following:
Fluiddynamic model
Microalgal cells are considered to be solid particles and fixed CO_{2} and released O_{2} are gases, therefore the analyzed system can be described as a multiphase flow and transport. However, it is possible to consider our multiphase system as a suspension by neglecting the gaseous phase. Then, the microalgal cells represent the dispersed phase of a suspension. Also, considering the cell density about 10 kg m ^{−3} for the dry weight of biomass, i.e., 1% of mass content, allows us to classify our system as a singlephase system based on Newtonian viscosity relationship within the employed computational software ANSYS Fluent [20]. Moreover, we assume that the mean diameter of a spherical microalgae cell is approximately ten micrometers [5] which supports this simplification.
Mass density of the suspension is determined as ρ=ρ_{w} (1−k)+ρ_{s} k, where ρ_{w} is the mass density of the medium and ρ_{s} represents the cell mass density, and k is the volume fraction. However, one can assume ρ=ρ_{w} because of a homogeneous distribution of algal cells (no aggregation) which are actually floating in the medium. Furthermore, the interparticle distances in our case of dilute suspension are sufficiently large to calculate flow field over each particle or cell [18], i.e., particles do not interfere with the flow field. This idea has been proposed by J. Pruvost et al. [21] who showed that the smallest eddy size, based on Kolmogorov scale, is approximately ten times bigger than the cell size, i.e., neglecting the cellvortex interaction.
Therefore, we model the incompressible liquid phase (suspension of water, nutrients and microalgae), thus the classical system of NavierStokes equations and the continuity equation is employed as:
in [ t_{0},T],×Ω, with suitable boundary conditions on [ t_{0},T],×∂Ω and initial conditions in Ω. v,p,f,ρ and ν denote the fluid velocity, the pressure, the body forces, fluid density and kinematic viscosity, respectively.
Reaction model and its timescale analysis
The photosynthetic reactions depend on many variables such as irradiance, accessibility of both inorganic (Fe, P, etc) and organic nutrients (e.g., carbohydrates), concentration of gases, osmolality, and many others. In this work, we focus only on the key variable determining the reaction kinetics, the fluctuations of irradiance; the other variables are assumed to be externally controlled to not limit the cellular growth. Generally, both light attenuation and scattering should be considered to describe the irradiance distribution inside MCS [16]. Moreover, for complicated MCS geometries and multiple light sources distribution, the irradiance field has to be determined experimentally. Here, for the sake of simplicity and manageable computational costs, we suppose the BeerLambert law (5) is valid, i.e., the incident irradiance I_{0} is exponentially decreasing with the increasing light path s:
where Λ is an attenuation coefficient describing the attenuation of irradiance by the suspension of microbial cells in the liquid medium (unit: m ^{−1}).
Based on the decades of experimental research of photosynthesis, the reliable model of reaction kinetics in MCS should implement no less than three time scales across the photosynthetic reactions, see [22]. The following systems are integrated in the reaction term R(c_{i}) in (1): (i) activation of the light harvesting complex (photosynthetic unit – PSU) in light reactions, (ii) biomass output, and (iii) photoinhibition, i.e., damage of pigmentprotein complexes by excessive irradiance. Photoacclimation (light adaptation) is another important process which, however, can be observed in the tens of hoursdays time scale. Therefore, it is possible to model it as a (quasi) steady state mechanism, e.g., implicitly via a reaction model parameter when another differential equation describes dynamics of this parameter [23].
Comparison of the characteristic times of microalgae growth (∼ hours) with the mixing due to dispersion – turbulent diffusion (∼ seconds) indicates a huge gap between both processes. There is a possibility to adopt a very simple growth modeling approach (e.g., Monod or Haldane kinetics) leading to two alternatives. The first one is governed by the time scale of chosen steady state kinetics and has to neglect details concerning mixing phenomena. The second alternative, based on the time step in order of seconds, can observe transitions because of the hydrodynamic mixing but has to neglect the changes in photosynthetic growth (in form of a growth kinetic model). However, both alternatives entirely miss the connection across the transport and reaction processes.
As a suitable candidate for the reaction model, we adopted a mechanistic threestate model of photosynthetic factory (PSF) [24–29]. PSF model considers 3 states in which microalgae cells may exist: activated – A, inhibited – B and rested – R. Although the fluiddynamical properties of cells in each of three states are exactly the same, these states are impacted differently by dynamic changes of the environment, i.e., spatiotemporal changes of states concentrations determine the biomass production, see (6).
Let be the concentrations of respective components c_{A}, c_{B}, and c_{R} (identical units as for the microalgae cell density c_{x} in whole PBR – generally 10^{6} cell ml ^{−1} as in [28]). Then the following relation holds (for ∀t∈[t_{0},t_{∞}], and ∀x∈∂Ω): c_{A}(x,t)+c_{B}(x,t)+c_{R}(x,t)=c_{x}(x,t). The dimensionless scalars y_{A}=c_{A}/c_{x}, y_{B}=c_{B}/c_{x} and y_{R}=c_{R}/c_{x} (molar fractions) are respective states of the PSF model, see the next section.
We note that the spatiotemporal averaged rate of photosynthetic production (specific growth rate \(\mu := \dot {c_{x}}/c_{x}\)) is related (constant κ) to the activated state fraction [28, 30]:
where term Me represents a cellular maintenance; Me is considered constant although it may vary if hydrodynamical shear stress is considered [29].
Decimal value 10^{−4} [s ^{−1}] of term κγ in (6) allows the shift from the time scale of light fluctuation (caused by fluid dynamics or the light source via the irradiance distribution I(x) in PSF model) to the time scale of biomass growth (macroscale) with no impact on the accuracy. Then, the integral in (6) is evaluated exactly and factor κγ in (6) provides the value of a specific growth rate. Thus, the employed PSF model satisfies the requirement for the time scales of the model candidate. Three of its five parameters can be seen as time constants for three different time scales: i) the light and dark reactions (∼ seconds), ii) the photoinhibition (∼ minutes) and iii) the microalgae growth (∼ hours). The remaining two parameters represent the optimal irradiance and the conversion from dimensionless growth to the real values, see [27, 31] for more details.
Model of photosynthetic factory – PSF model
Here we characterize the multi timescale threestate model of photosynthetic factory (PSF model) proposed by Eilers and Peeters [24, 25] and further developed by Wu and Merchuk [28, 29] and Celikovsky, Papacek and Rehak [26, 27, 30–32]. PSF model (Fig. 1), is used for derivation of the reaction term R(c_{i}), see the transport Eq. (1).
PSF model (7) incorporates the dynamics of three fundamental states propagating in different orders of time scales: (i) cell growth (including the shear stress effect), (ii) photoinhibition, and (iii) photosynthetic light and dark reactions. The state vector of the PSF model is three dimensional, y=(y_{R},y_{A},y_{B})^{⊤}, where y_{R} indicates the probability that PSF is in the resting state R, y_{A} the probability that PSF is in the activated state A, and y_{B} the probability that PSF is in the inhibited state B.
The values of PSF model parameters, i.e., α, β, γ, δ, see (7), and κ, see (6), are taken from original study [28], where Wu and Merchuk identified these values of PSF model parameters for the microalga Porphyridium sp.: α = 1.935 × 10^{−3} μE^{−1} m^{2}, β = 5.785 × 10^{−7} μE^{−1} m^{2}, γ = 1.460 × 10^{−1} s ^{−1}, δ = 4.796 × 10 ^{−4} s ^{−1}, κ = 3.647 10 ^{−3} and Me = 0.059 h ^{−1}. For details regarding the experimental setup, for parameter estimation and the identifiability study, see [31, 32].
In order to make the PSF model parameter estimation more robust, the following reparametrization was introduced and the singular perturbation method was used in [31]: \(q_{1}: = \sqrt {\frac { \gamma \delta }{ \alpha \beta }}, ~ q_{2}: = \sqrt {\frac {\alpha \beta \gamma }{\delta {(\alpha +\beta)}^{2}}},~ q_{3} := \kappa \gamma \sqrt {\frac { \alpha \delta }{\beta \gamma }},~ q_{4}:=\alpha q_{1},~ q_{5} := \beta /\alpha.\) Consequently, the PSF model acquires the following form:
For the given input variable, i.e., the irradiance u(t), the ODE system (7) can be solved either by numerical methods or by asymptotic methods. However, the irradiance level of 250 μE m ^{−2} s ^{−1} (the value of the parameter q_{1} which maximizes the steadystate growth [31]), provides two negative eigenvalues of the (7) system matrix: λ_{1} = 0.63, λ_{2} = 0.59 10 ^{−3}. The ODE system (7) is stiff due to value 10^{3} of the ratio \(\frac {\lambda _{1}}{\lambda _{2}}\). This fact points to the existence of two processes: (i) photosynthetic reactions, and (ii) photoinhibition, widely separated from the point of view of their characteristic times. Further, without significant loss of accuracy, we use the socalled fast reduction [27], i.e., the behavior of the system (7) is characterized by the only one following ODE
and the “slow” variable y_{B} can be regarded as a constant depending on the averaged value u=u_{av}. According to our works [27, 33] it holds
Both the nonreduced (7) and the reduced order PSF model Eq. (11) have been incorporated as a UserDefined Function (UDF) in ANSYS Fluent, see the following section.
Numerical aspects of PDEs implementation and solution
The key issue is in the incorporation of a CFD code and photosynthetic reaction kinetics in one modeling framework. A proper analysis of characteristic times of microalgae growth and of mixing due to the convective motion is especially important. Two appealing approaches, Lattice Boltzmann method [34] and hybrid multicompartment/CFD approach [35], were previously used for bioreactor modeling [36, 37]. Both methods seem to be computationally efficient, however, the discussion of their advantages and disadvantages is far of the scope of this paper. Our approach is based on the implementation of a UserDefined Function (UDF) within the commercial CFD code ANSYS Fluent, which provides the possibility to define an arbitrary reaction term, see below.
As denoted in the previous section, ODE system (7) can be simplified to only one differential Eq. (11). The righthand part of this equation indicates the rate of change of activated state y_{A} which can be used in the definition of UDF function in ANSYS Fluent. Macro DEFINE_VR_RATE provides the possibility to define the arbitrary reaction term in (1), including its dependency on the spatial coordinate as represented by (5). Macro C_CENTROID can retrieve the corresponding coordinates of the current mesh cell. Three states of the microalgae culture were represented as individual species with the same molar weights in ANSYS Fluent. Therefore, mass and molar fractions should be identical. Example of the userdefined reaction rate definition employed in our simulations for the reduced system (11) follows:
We performed a grid independence study for three mesh sizes of the given geometry (2D square cavity with moving wall), see Fig. 2. Grid Convergence Index (GCI) describing an average velocity of the flow field for the finest mesh was evaluated as 1.14% according to Celik et al. [38]. All simulations of the 2D geometry were performed on a grid with approx. 10 thousand mesh elements. Refined mesh elements near walls were created so that dimensionless thickness Y^{+} was around 1 to capture properly the viscous sublayer in the turbulent flow regime. Time step 0.025 s was used and the accuracy of the solution at the end of time range 0 – 5 s was evaluated as 0.19%.
Every simulation consisted of two steps. In the first step, the NavierStokes Eq. (4) were solved iteratively to get a converged solution of the steadystate flow field for the corresponding mixing rate (Reynolds number). In the second step, the flow field was assumed to be fully developed and only Eq. (1) describing the transport and reaction of individual species (states of the algae in our case) were solved. This approach is more efficient but in cases when the flow field exhibits a transient behaviour, the NavierStokes equations must be solved simultaneously with the species transport equations. For example, the CoutteTaylor device shows an inherently transient behaviour at some specific rotations rates [39]. Periodic oscillations can be observed even in the square cavity [40, 41] which we neglected in this study.
Results and discussion
Our early works [27, 30] theoretically confirmed the phenomenon of flashing light enhancement, known from experiments [13–15]. However, these models were applicable only for specific PBR geometries, being a CouetteTaylor bioreactor in [30] and flat panel PBR in [27], and suffered from oversimplification (NavierStokes equations were not solved there within a fluid dynamic model). Here, for the first time, we introduce a general approach applicable to any geometry of MCS and thus ready for industrial applications. For the validation of our integrated model, we have chosen the special geometry consisting of a 2D square cavity with one moving wall (Fig. 2), representing an idealized part of a CouetteTaylor device, see Fig. 3. We note that microalgae cells are forced to travel between dark and light sides of the CouetteTaylor device (see Taylor vortex in [39]) which is an analogy to the flashing light regime and consequently the flashing light enhancement [15]. Indeed, very high cell densities (about 150 g/l) were reported in this type of PBR by E. A. Davis in [42]. The author claims that the growth enhancement is due to an ordered mixing (E. A. Davis misuses the term “turbulence”) produced by the CouetteTaylor device. Nevertheless, the mixing enhanced light and mass transfer is only one part of the whole picture. Based on the fact that the standing issue of microalgae mass production in PBR is the cell fragility [43], we see also the other enhancing effect  low shear stress transmitted to microbial cells in the laminar Taylor vortex flow regime [44].
Photosynthesis in fluctuating light regime
The growth of photosynthetic microorganisms could be described by the steadystate photosynthetic reaction kinetics called P–I (photosynthesisirradiance) curve representing the photosynthetic response in small cultivation systems with a homogeneous light distribution [17]. However, MCS models based on P–I curves only do not take into account the effect of hydrodynamic mixing (or mixing via static mixers) which has a different characteristic time scale compared to the time scale of the steadystate reaction kinetics, see “Reaction model and its timescale analysis” section and [27, 31] for more details. In order to overcome this problem, Terry in [15] employed “artificial” corrections to the steadystate kinetic model reflecting the enhancement in the photosynthetic efficiency during fluctuating light conditions present in the dynamic model.
In our work, we aim to build a general predictive model, independent of specific experimental data. Such model must be built on universal principle that timescales of all relevant processes are considered; see details in section “Reaction model and its timescale analysis”.
To determine if our modeling framework is able to describe the flashing light enhancement phenomenon, we have chosen a 2D square cavity with the moving top wall as our hypothetical MCS, see Fig. 2. This MCS can be viewed as the representation of a CouetteTaylor device, i.e., one cell from its axial crosssection when the socalled Taylor vortex flow regime takes place. This occurs when the Taylor number exceeds the first critical value [39], see Fig. 4. The cultivation system is illuminated from the bottom side only by the incident irradiance I_{0} and the irradiance level is exponentially decreasing with the increasing zcoordinate until the top (z=L) according to the BeerLambert law (5). It is convenient to work further with the normalized irradiance u=I/I_{opt}, where I_{opt} reaches the value of 250 μE m ^{−2} s ^{−1} [28]. Our analysis showed that the value of the averaged irradiance within MCS u_{av}=1, see Table 1, provides an optimal condition for the microalgae culture growth. This optimization problem depends mainly on mixing intensity and it is discussed elsewhere [33]. Only for the lumped parameter system, the value u=1 is the optimal one. The relations for the irradiance level in depth z and for the average (absorbed) irradiance in the whole MCS are then:
where u_{0}=I_{0}/I_{opt} is the incident irradiance on the bottom (z=0) and L is the MCS thickness in direction of light gradient. It is further convenient to define a dimensionless “optical thickness” of the system Ot>0 as follows: Ot:=ΛL. It holds \(\frac {r_{1/2}}{L}=\frac {\ln (2)}{Ot}\), where r_{1/2} is the length interval (unit: m) in which the light intensity diminishes to one half. Then, ratio \(\frac {u(z)}{u_{av}}\) for optically thin culture system (Ot≈5.5 in our case) has
We suggest that the “irradiance history” of an individual cell within our system in Fig. 2 corresponds to the hydrodynamically induced light/dark cycles shown in Fig. 3 (produced by the CouetteTaylor device). The irradiance history can be calculated by linking the trajectory of an individual microalgae cell and the irradiance field within device.
The final step is the implementation of the transportreaction model based on the PDEs (14). The key role plays the reaction term R(c_{i}) in (1). The R(c_{i}) term is “local” in sense that it depends on the actual transport and reaction rates, and it contains all relevant time scales, thus, it perfectly fits to our requirements. The technical simplicity of PDE based approach, which takes advantage of the sophisticated CFD codes, e.g., ANSYS Fluent [20], is obvious when compared to the alternative Lagrangian approach. While the Eulerian approach reaches its governing equations from the mass balance in an infinitesimal control volume, the Lagrangian approach is based on the above mentioned irradiance history of an individual cell. Although conceptually simple, the application of Lagrangian approach to MCS modeling presents serious problems. First, either an experimental technique, e.g., computerautomated radioactive particle tracking (CARPT) [45] or a CFD code is needed for the cell trajectory description, see Fig. 3. Second, using an ensemble of cell trajectories across the irradiance field inside MCS, the probabilistic description of “irradiance history” of individual cells is reached. Finally, the problem reduces to a reaction kinetics model represented by an ordinary differential equation with a stochastic input.
The employed Eulerian approach is more convenient, mainly because it benefits from CFD codes and the implementation of an appropriate dynamic model of microbial growth (with lumped parameters) into a distributed parameter system is merely technical.
Case study: Biological performance of MCS in flashinglight regime
Having the modeling framework assembled, we aimed to reproduce the hydrodynamic regime where high frequency lightdark cycles are induced and trigger the flashing light enhancement [14, 15]. The values of parameters employed in the simulation of MCS performance are summarized in Table 1.
The PSF model consists of the ODE system (7), nevertheless it can be reduced by using the singular perturbation method [33]. Then, the behavior of the system (7) is characterized by only one ODE (11). Furthermore, for the local description of PSF model, we need the expression for u(z) in the case of average irradiance u_{av}=1; based on (14): \(u(z)\approx \Lambda L~{e^{\Lambda L \frac {z}{L}}}\).
For the initial condition and other parameters in Table 1, the case study was resolved using CFD code ANSYS Fluent: PDE (14) embedded and the reaction kinetics implemented as a special UDF (UserDefined n) according to the description in section “Numerical aspects of PDEs implementation and solution”. The resulting steadystate solution for molar fraction y_{A} is shown in Fig. 5. The spatial distribution of molar fraction y_{B} steadystate solution is shown in Fig. 6. The spatiotemporal average for molar fraction y_{A} is shown in Fig. 7 and, finally, the temporal evolution of y_{B} steadystate is depicted in Fig. 8 (right). Notice nice patterns of y_{B} for low mixing conditions in Fig. 6 (left and center).
It is clear that the presented simulations are in the qualitative match with the previously published experimental data [14], see Fig. 7 on the left. Also, the time scale of y_{A} state dynamics is shown in Fig. 8, describing the time course of the spatial average of y_{A} for three different velocities v_{L}, corresponding to Reynolds number of 0, 1000, and 100000, respectively. The mixing intensity determines the frequency of hydrodynamically induced light/dark cycles. Similarly to our previous simplified studies [27, 30], we can see the hydrodynamically induced growth enhancement due to the flashing light. The major difference and improvement to previous studies is that presented modeling framework can be used to any MCS setup and not only to special cases.
Analysis of our numerical model indicates that the process is characterized by the Damköhler number Da:=t_{tr}/t_{r} which describes the ratio between the characteristic time of transport (\(t_{tr}=\frac {L}{v_{L}}\)) and reaction time (t_{r}=1/[q_{4}(u_{av}+q_{2})]). In chemical engineering, the Damköhler numbers Da are frequently used to describe the ratio between the chemical reaction time scale (reaction rate) and the flow time scale (convective mass transport rate) occurring in a system. For our specific case, when the characteristic time t_{r} is fixed, it holds:
It is known that extremely high or low values of the characteristic numbers leads to numerical difficulties during solving PDEs, e.g., (1). Employing Monod or Haldane kinetics (steadystate) might lead to such difficulties and that is why a multi timescale model would be more appropriate.
In order to quantify the impact of mixing velocity on the cellular growth, we defined an objective function J_{MCS} as a volumetric productivity (usually in grams per liter per day):
For a usual case of steady state operation, it is only the activated state fraction y_{A}(u,z,t) in J_{MCS} which has to be integrated over the MCS domain. Thus, in Fig. 7 we show the dependence of the normalized performance index
We note that the empirical data have an illustrative and testing purpose only and for quantitative simulations one has to employ parameters for specific MCS and particular microalgae. However, the result of biological performance for our MCS compared to experimental results shown in Fig. 7 (left) clearly shows that the positive impact of mixing on the growth induced by the convective motion, i.e., the simulated flashing light enhancement is comparable to experiments. Analyzing simulations in Fig. 7 (right) implies that lower Da (bigger v, better mixing) leads to better performance J (approaching its maximal theoretical value 0.625). Nevertheless, the harmful impact of the hydrodynamically induced shear stress originally considered in Me term [29], see (6), is not taken into account and thus in reality certain mixing intensity would start damaging the cellular integrity.
Conclusions
In order to reach commercially sustainable production for algae biofuels, the empirical data indicates the need for accurate and reliable modeling approach with both predictive and optimizing capabilities. In this work, we presented the general unified modeling framework for microalgae culture systems (MCS) which can be applied to an arbitrary geometry of the production system and any microbial strain in general. This allows evaluating the performance of various MCS architectures, operating conditions, and microalgae strains as well. All parts of the multidisciplinary framework, i.e., the growthrate model, the fluiddynamic model, and the irradiancedistribution model, are integrated within CFD package ANSYS Fluent. In the illustrative case study, we have shown that a threestate model of photosynthetic factory (PSF model) well behaves under hydrodynamically induced high frequency lightdark cycles regime and copes with the requirement imposed on the reaction model, i.e., it correctly describes both the quasi steadystate (flashing light enhancement) and dynamic phenomena (time course of the PSF model activated state).
We can conclude that our modeling framework provides an adequate and physically accurate description of microalgae growth in a MCS and thus it is suited for the optimal control problem formulation as well. To the extent of our knowledge, this is the first time when such advanced solution was successfully used in biotechnology of microalgae.
Our future plans are to further focus on the optimization problem and its formulation. Finally, we plan to analyze the shear stress impact on the microalgae growth.
Abbreviations
 CARPT:

Computerautomated radioactive particle tracking
 CFD:

Computational fluid dynamics
 GCI:

Grid convergence index
 MCS:

Microalgae culture systems
 PBR:

Photobioreactors
 PDEs:

Partial differential equations
 PI:

Photosynthesisirradiance
 PSF:

Photosynthetic factory
 UDF:

Userdefined function
References
 1
Su Y, Song K, Zhang P, Su Y, Cheng J, Chen X. Progress of microalgae biofuel’s commercialization. Renew Sust Energ Rev. 2017; 74:402–11. https://doi.org/10.1016/j.rser.2016.12.078.
 2
Burlew J. Algal Culture From Laboratory to Pilot Plant. Washington: Carnegie Institution of Washington; 1953.
 3
Grobbelaar JU. Microalgal biomass production: challenges and realities. Photosynth Res. 2010; 106(1):135–44. https://doi.org/10.1007/s1112001095735.
 4
Ooms MD, Dinh CT, Sargent EH, Sinton D. Photon management for augmented photosynthesis. Nat Commun. 2016; 7:12699. https://doi.org/10.1038/ncomms12699.
 5
Richmond A. Handbook of Microalgal Culture: Biotechnology and Applied Phycology. Oxford: Blackwell Publishing; 2008.
 6
Chisti Y. Biodiesel from microalgae, biotechnology advancesbiodiesel from microalgae. Biotechnol Adv. 2007; 25:294–3062.
 7
Bernardi A, Perin G, Sforza E, Galvanin F, Morosinotto T, Bezzo F. An identifiable state model to describe light intensity influence on microalgae growth. Ind Eng Chem Res. 2014; 53(16):6738–49. https://doi.org/10.1021/ie500523z. PMID: 25678739. https://doi.org/10.1021/ie500523z.
 8
Park S, Li Y. Integration of biological kinetics and computational fluid dynamics to model the growth of nannochloropsis salina in an open channel raceway. Biotech Bioeng. 2015; 112(5):923–33. https://doi.org/10.1002/bit.25509.
 9
Sato T, Yamada D, Hirabayashi S. Development of virtual photobioreactor for microalgae culture considering turbulent flow and flashing light effect. Energy Convers Manag. 2010; 51:31196–201.
 10
MullerFeuga A, Guédes RL, Pruvost J. Benefits and limitations of modeling for optimization of Porphyridium cruentum cultures in an annular photobioreactor. J Biotechnol. 2003; 103:153–63.
 11
Bitog JP, Lee IB, Lee CG, Kim KS, Hwang HS, Hong SW, Seo IH, Kwon KS, Mostafa E. Application of computational fluid dynamics for modeling and designing photobioreactors for microalgae production: A review. Comput Electron Agric. 2011; 76(2):131–47. https://doi.org/10.1016/j.compag.2011.01.015.
 12
Bernard O, Mairet F, Chachuat B In: Posten C, Feng Chen S, editors. Modelling of Microalgae Culture Systems with Applications to Control and Optimization. Cham: Springer: 2016. p. 59–87. https://doi.org/10.1007/10_2014_287.
 13
Kok B. Experiments on photosynthesis by chlorella in flashing light In: Burlew J, editor. Algal Culture From Laboratory to Pilot Plant, vol. 600. Washington: Carnegie Institution of Washington: 1953. p. 63–75.
 14
Nedbal L, Tichý V, Xiong F, Grobbelaar JU. Microscopic green algae and cyanobacteria in highfrequency intermittent light. J Appl Phycol. 1996; 8:325–33.
 15
Terry KL. Photosynthesis in modulated light: Quantitative dependence of photosynthetic enhancement on flashing rate. Biotech Bioeng. 1986; 28:988–95.
 16
Cornet JF, Dussap CG, Gros JB, Binois C, Lasseur C. A simplified monodimensional approach for modeling coupling between radiant light transfer and growth kinetics in photobioreactors. Chem Eng Sci. 1995; 50:1489–500.
 17
Schugerl K, Bellgardt KH. Bioreaction Engineering, Modeling and Control. Berlin Heidelberg: Springer; 2000.
 18
Beek WJ, Muttzal KMK, Heuven JWV. Transport Phenomena, 2nd edn. Chichester: Wiley; 1999.
 19
AlvarezV ázquez LJ, Fernández FJ. Optimal control of bioreactor. Appl Math Comput. 2010; 216:2559–75.
 20
Fluent ANSYS. Theory guide. Ansys Inc. 2015.
 21
Pruvost J, Legrand J, Legentilhomme P, MullerFeuga A. Simulation of microalgae growth in limiting light conditions: Flow effect. AIChE J. 2002; 48(5):1109–20. https://doi.org/10.1002/aic.690480520.
 22
Han BP, Virtanen M, Koponen J, Straškraba M. Effect of photoinhibition on algal photosynthesis: a dynamic model. J Plankton Res. 2000; 22(5):865. https://doi.org/10.1093/plankt/22.5.865.
 23
Kmeť T, Straškraba M, Mauersberger P. A mechanistic model of the adaptation of phytoplankton photosynthesis. Bull Math Biol. 1993; 55:259–75.
 24
Eilers PHC, Peeters JCH. A model for the relationship between light intensity and the rate of photosynthesis in phytoplankton. Ecol Model. 1988; 42(3):199–215. https://doi.org/10.1016/03043800(88)900579.
 25
Eilers P, Peeters J. Dynamic behaviour of a model for photosynthesis and photoinhibition. Ecol Model. 1993; 69(12):113–33.
 26
Papacek S, Celikovsky S, Stys D, RuizLeón J. Bilinear system as modelling framework for analysis of microalgal growth. Kybernetika. 2007; 43:1–20.
 27
Papacek S, Matonoha C, Stumbauer V, Stys D. Modelling and simulation of photosynthetic microorganism growth: random walk vs. finite difference method. Math Comput Simul. 2012; 82(10):2022–32. https://doi.org/10.1016/j.matcom.2011.07.006. The Fourth IMACS Conference : Mathematical Modelling and Computational Methods in Applied Sciences and Engineering, Devoted to Owe Axelsson in ocassion of his 75th birthday.
 28
Wu X, Merchuk JC. A model integrating fluid dynamics in photosynthesis and photoinhibition processes. Chem Eng Sci. 2001; 56:3527–38.
 29
Wu X, Merchuk JC. Simulation of algae growth in a benchscale bubble column reactor. Biotech Bioeng. 2002; 80(2):156–68. https://doi.org/10.1002/bit.10350.
 30
Papacek S, Stumbauer V, Stys D, Petera K, Matonoha C. Growth impact of hydrodynamic dispersion in a couette–taylor bioreactor. Math Comput Model. 2011; 54(7):1791–5. https://doi.org/10.1016/j.mcm.2010.12.022. Mathematical models of addictive behaviour, medicine & engineering.
 31
Rehak B, Celikovsky S, Papacek S. Model for photosynthesis and photoinhibition: Parameter identification based on the harmonic irradiation O _{2} response measurement. TAC IEEE and TCAS IEEE. 2008:101–108.
 32
Papacek S, Celikovsky S, Rehak B, Stys D. Experimental design for parameter estimation of two timescale model of photosynthesis and photoinhibition in microalgae. Math Comput Simul. 2010; 80(6):1302–9. https://doi.org/10.1016/j.matcom.2009.06.033. Fifth IMACS Seminar on Monte Carlo Methods Applications of Computer Algebra 2007 (ACA 2007) special session on Nonstandard Applications of Computer Algebra Computational Biomechanics and Biology, a collection of papers presented at the 1st IMACS International Conference on the Computational Biomechanics and Biology ICCBB 2007.
 33
Celikovsky S, Papacek S, CervantesHerrera A, RuizLeon J. Singular perturbation based solution to optimal microalgal growth problem and its infinite time horizon analysis. IEEE Trans Autom Control. 2010; 55(3):767–72. https://doi.org/10.1109/TAC.2010.2040498.
 34
Succi S. The Lattice Boltzmann Equation. New York: Oxford University Press; 2001.
 35
Bezzo F, Macchietto S, Pantelides CC. Computational issues in hybrid multizonal/computational fluid dynamics models. AIChE J. 2005; 51(4):1169–77. https://doi.org/10.1002/aic.10383.
 36
Bezzo F, Macchietto S, Pantelides CC. General hybrid multizonal/cfd approach for bioreactor modeling. AIChE J. 2003; 49(8):2133–48. https://doi.org/10.1002/aic.690490821.
 37
Štumbauer V, Petera K, Štys D. The lattice Boltzmann method in bioreactor design and simulation. Math Comput Model. 2013; 57:1913–8.
 38
Čelikovský S. Procedure for estimation and reporting of uncertainty due to discretization in cfd applications. J Fluids Eng. 2008; 130(7). https://doi.org/10.1115/1.2960953.
 39
Taylor GI. Stability of a viscous liquid containing between two rotating cylinders. Phil Trans R Soc. 1923; A223:289–343.
 40
Goodrich JW. An unsteady timeasymptotic flow in the square driven cavity. Technical report tech. mem. 103141. NASA. 1990.
 41
Bruneau CH, Saad M. The 2d liddriven cavity problem revisited. Comput Fluids. 2006; 35:326–48.
 42
Davis EA. Turbulence In: Burlew J, editor. Algal Culture From Laboratory to Pilot Plant. Washington: Carnegie Institution of Washington: 1953. p. 135–138.
 43
Gudin C, Chaumont D. Cell fragility — the key problem of microalgae mass production in closed photobioreactors. Bioresour Technol. 1991; 38(2):145–51. https://doi.org/10.1016/09608524(91)90146B.
 44
Haut B, Amor HB, Coulon L, Jacquet A, Halloin V. Hydrodynamics and mass transfer in a couette–taylor bioreactor for the culture of animal cells. Chem Eng Sci. 2003; 58(3):777–84. https://doi.org/10.1016/S00092509(02)006073. 17th International Symposium of Chemical Reaction Engineering (IS CRE 17).
 45
Luo HP, Kemoun A, AlDahhan MH, Sevilla JMF, Sánchez JLG, Camacho FG, Grima EM. Analysis of photobioreactors for culturing highvalue microalgae and cyanobacteria via an advanced diagnostic technique: Carpt. Chem Eng Sci. 2003; 58(12):2519–27. https://doi.org/10.1016/S00092509(03)000988.
Funding
This work was supported by the Ministry of Education, Youth and Sports of the Czech Republic  projects ”CENAKVA” (No. CZ.1.05/2.1.00/01.0024), ”CENAKVA II” (No. LO1205 under the NPU I program) and ’The CENAKVA Centre Development’ (No. CZ.1.05/2.1.00/19.0380). Publication costs were funded by University of South Bohemia.
About this supplement
This article has been published as part of BMC Systems Biology Volume 12 Supplement 5, 2018: Selected articles from the 5th International WorkConference on Bioinformatics and Biomedical Engineering: systems biology. The full contents of the supplement are available online at https://bmcsystbiol.biomedcentral.com/articles/supplements/volume12supplement5.
Author information
Affiliations
Contributions
SP designed the model and drafted the manuscript, KP performed the simulations and collaborated in drafting the manuscript. JJ validated the results and supervised the manuscript preparation. All authors interpreted part of the data and approved the final manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
Papacek, S., Jablonsky, J. & Petera, K. Advanced integration of fluid dynamics and photosynthetic reaction kinetics for microalgae culture systems. BMC Syst Biol 12, 93 (2018). https://doi.org/10.1186/s1291801806119
Published:
Keywords
 Microalgae
 Mathematical modeling
 Photosynthesis
 CFD
 Microalgae culture systems
 Flashing light enhancement