A proof of Kirchhoff's first law for hyperbolic conservation laws on networks

1.   Introduction

During the last few years, third generation biorefinery processes, converting waste materials into multiple valuable products, have attracted increasing attention. In particular, industries and researches are focusing on integrated biological systems for the production of bioenergy and biopolymers. This new waste management approach is environmental-friendly and makes two contributions toward reducing greenhouse gas (GHG) emissions and replacing fossil fuels ^[1].

Among different processes and configurations, photo fermentation (PF) represents one of the most suitable alternatives. Indeed, Purple Non Sulfur Bacteria (PNSB) are able to contextually produce hydrogen (H₂) and polyhydroxybutyrate (PHB) in a single stage process, therefore reducing production costs and management issues ^[2]. In addition, PF provides almost complete degradation of organic compounds from the treated wastewater ^[3]. However, PF is affected by many process conditions (e.g., type and concentration of carbon and nitrogen source ^[4,5,6], light availability ^[4,6,7], temperature ^[6,7], pH ^[6,7], presence of inhibiting substances ^[5,6]) which need to be considered for the effective optimization of hydrogen and PHB production. Moreover, more studies are still required to completely clarify the process mechanism and its efficiency under different environmental conditions.

In this context, the use of mathematical models is helpful to effectively control the PF process. Mathematical models can simulate the influence of different environmental and operational conditions and limit the load of the experimental tests, which are costly and time demanding ^[8,9]. In addition, the application of mathematical modeling to complex biological processes, such as fermentation ^[10], microbial electrolysis ^[11], anaerobic digestion, and bio photolysis ^[12], allows to identify and explore the most relevant dynamical interactions of bacteria competing for substrates and space in the same reaction environment. These aspects are crucial for the correct management and/or design of bioreactors aimed at producing value-added chemicals or biogas ^[13,14].

Due to the multiple effects of the process conditions on biological kinetics, PF is a complex process to model ^[15]. Over the past few years, different methodologies have been adopted for the mathematical definition of the PF process. An extensive review of such models is provided in Policastro et al. ^[16]. The most frequently applied equations are based on the Monod's equation and its modification ^[17]. The Luedeking-Piret model, the Gompertz equation and their modifications have been adopted as well ^[18,19]. Recently, nitrogen (N) and phosphorous (P) recovery from domestic wastewater by PNSB has been described using a novel ADM1-based mechanistic model ^[20].

Despite the increasing attention of the scientific community on valuable biological routes to accumulate PHB and produce biohydrogen, there is a lack of mathematical models describing the concomitant hydrogen and PHB generation achievable in PF processes. Therefore, the aim of this work is the development of a novel and complete mathematical model able to describe the photofermentative bioconversion of organic compounds in both hydrogen and PHB, by PNSB. In particular, the model has been specialized for organic acids conversion by Rhodopseudomonas palustris, operating in batch conditions, and under continuous illumination conditions.

The model definition and its study were carried out based on the following steps:

1) based on biological evidence revealing PF mechanisms in natural and industrial environments, a system of ordinary differential equations able to reproduce the PF process has been defined;

2) a local sensitivity analysis for selected parameters has been performed as a supporting tool for the calibration of the mathematical model;

3) the model has been calibrated and validated using different lab scale experimental data sets;

4) the model predictions have been further investigated quantitatively to test the model consistency in describing relevant mechanisms;

5) finally, a fast-slow analysis technique has been carried out to investigate the additional H₂ production from accumulated PHB in bacterial cells.

2.   Photo fermentation mechanisms and process conditions

The model has been defined based on the current knowledge of the PF biological mechanisms. PNSB are able to concomitantly produce hydrogen and PHB when photoheterotrophic growth occurs. Under this condition, PNSB get electrons and carbon from reduced organic compounds, such as acetate and other organic acids, alcohols, and carbohydrates. The substrate is converted to hydrogen, other catabolites, new microorganisms and PHB intercellular granules ^[21]. The pathway toward hydrogen production is supported by the ATP-dependent nitrogenase enzyme. The required ATP is synthetized by using light as energy source. Therefore, light availability and intensity notably influence the PF process. Previous studies have showed that high light intensities result in the partial inhibition of the hydrogen production ability ^[22]. In addition, excessively high biomass concentrations avoid the penetration of light into the cultivation system due to the self-shading phenomenon ^[23]. Consequently, the decrease in light availability inhibits the ATP formation and reduces the hydrogen yield. The ATP formation can be inhibited by excessively alkaline or acidic pH ^[24]. The incubation temperature is another relevant factor, which influences bacteria growth, hydrogen production and PHB accumulation ^[21].

Substrate characteristics and concentration also has a major effect on PNSB activities. In particular, high substrate concentrations lead to the accumulation of acids in the reaction environment, inhibiting the process ^[25]. Moreover, the presence of high ammonia concentrations or toxic compounds (e.g., phenols) in the culture system decreases the nitrogenase activity leading to the partial or complete inhibition of PNSB ^[26]. The availability of nutrients is of primary concerns as well. In particular, nitrogen is required for cell synthesis. It has been demonstrated that, nitrogen availability can regulate growth phases of the PNSB specie Rhodopseumonas palustris ^[27]. In presence of both carbon and nitrogen sources, bacteria can perform high growth rates (exponential phase). Subsequently, the stationary phase occurs without significant nitrogen depletion, and the biomass concentration is maintained by the presence and consumption of the carbon sources ^[27]. Experimental evidence demonstrates that hydrogen production takes place in all growth phases. However, Rhodopseudomonas palustris use more energy for anabolic reactions during the exponential phase, while shift their metabolism to catabolic reactions during the stationary phase ^[27,28].

Concerning mechanisms leading to the contextual hydrogen and PHB production, usually PHB synthesis occurs along with hydrogen production, when microorganisms are forced to live under stress conditions, such as nutrient scarcity or other nonoptimal environmental conditions ^[29]. In Rhodopseudomonas palustris and mixed PNSB cultures, it has been observed that the PHB synthesis prevented bacteria from experiencing the stress condition, and retained H₂ production ^[30,31]. Moreover, PNSB use PHB as a reserve for their survival when the carbon substrate became scarce. Therefore, the PHB depletion is usually observed when organic substrates reach low concentration levels ^[31].

3.   Model assumptions and description

The PF process is governed by many biotic and abiotic conditions and the metabolic mechanisms occurring in photobioreactors are not completely elucidated. The definition of a mechanistic mathematical model able to adequately reproduce real PF data requires to adopt simplifications to avoid the over-parametrization of the problem. Several authors have developed different mathematical models, considering the effect of a limited number of factors ^[16]. The selection of these factors was usually based on the specific application. For instance, the effect of light intensity has been included by Gadhamshetty at al. ^[19], in the logistic equation and the Luedeking-Piret model. Moreover, Sevinc et al. ^[18] considered the light intensity effect, modifying the previously introduced first order kinetic. The self-shading phenomenon and the substrate inhibition have been taken into account in several previous models, using various kinetic equations ^[19,32]. However, other authors considered the effect of pH and temperature in Monod type kinetics ^[15]. In addition, Puyol et al. ^[20] considered the ammonia inhibition phenomenon in a mechanistic model describing the nutrient removal process from wastewater by PNSB.

In the present work, a mathematical model able to describe hydrogen and PHB production by Rhodopseudomonas palustris at a lab scale is presented. The model is formulated for a stirred tank reactor fed in batch mode to reproduce the most common experimental set-up of PF literature. The concentration of this PNSB species has been related to the contextual production of H₂ and PHB by using a classical Monod-Michaelis-Menten scheme regulated by substrate availability. The general ODEs framework adopted for the evolution of a generical substrate ($S$), microbial species ($X$), and product ($P$) over time is written as:

where ${\nu }_{max}$ is the degradation rate of the substrate $S$, ${K}_{S}$ is the semi saturation constant, and $Y$ represents the yield of the biomass on the substrate.

The model included the most relevant factors influencing the biological process during PF batch tests, such as the self-shading and the substrate inhibition phenomenon. The ammonia inhibition and the presence of toxic compounds influencing Rhodopseudomonas palustris metabolism have been neglected. Due to the operating conditions usually adopted for PF indoor batch tests, a constant pH, temperature, and light intensity level have been considered. Consequently, the effects of these factors on PF performance have been omitted. To allow the model calibration and the chemical oxygen demand (COD) balance closure, the variable "catabolites" has been added to the model. This variable accounts for all the residual catabolites which were not detected by experimental analysis. The Rhodopseudomonas palustris growth has been described by using an original inhibition-like function able to regulate the biomass yield when the nitrogen source is relatively low. The model takes into account the concomitant PHB accumulation and hydrogen production due to microbial activity. Under substrate starvation, the model is able to account for the typical PHB depletion occurring in PF experiments. Indeed, the consumed PHB can be used for biomass maintenance and to produce an additional fraction of hydrogen and other catabolites.

4.   Model equations and numerical integration

The mathematical model consists of a system of 6 ordinary differential equations (ODEs). The equations describe the biomass growth, the substrates consumption, the PHB accumulation, the hydrogen production, the production of other catabolites, and the nitrogen depletion connected to the biomass growth. Nonlinear Monod-like kinetic functions have been introduced, and a first order kinetic has been used for biomass decay. The biomass growth and the substrate and nitrogen depletion are described by Eqs (4)–(6):

where $X$, $N$, and ${S}_{i}, i = 1, \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }2, \mathrm{ }3,$ are the microbial biomass, the nitrogen, and the substrates concentration, respectively. $Y$ denotes the biomass yield and $I$ are limiting/inhibition functions. ${f}_{i}$ are the Monod type kinetics, where $i = 1, \mathrm{ }2, \mathrm{ }3$ represent 3 different organic substrates for PNSB. ${f}_{4}$ is the Monod-like function for PHB concentration. Finally, ${K}_{d}$ is the biomass decay coefficient, and ${N}_{bac}$ is the nitrogen content of bacteria. Eqs (4)–(6) describe the PHB evolution, the hydrogen production, and catabolites accumulation over time as:

where ${F}_{i}$ denotes the hydrogen yield on the different substrates, ${F}_{{cat}_{i}}$ denotes the hydrogen yield on catabolites, and ${F}_{p}$ represents the hydrogen yield on PHB. Note that Eq (7) is defined by two different terms: the first describes the PHB consumption, and the second one accounts for the PHB accumulation in PNSB cells. The Monod functions ${f}_{i}, i = 1, \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }2, \mathrm{ }3$, and ${f}_{4}$ have been expressed as:

where ${\nu }_{{i}_{max}}$ and ${\nu }_{{P}_{max}}$ are the substrates and PHB consumption rates, respectively. ${K}_{{S}_{i}}, i = 1, \mathrm{ }\mathrm{ }\mathrm{ }\mathrm{ }2, \mathrm{ }3$, and ${K}_{{P}_{1}}$ are the half saturation constants for substrates and PHB, respectively.

As reported in the ADM1 (Batstone et al., 2022), non-competitive inhibition functions have been applied for inhibitory phenomena. Specifically, ${I}_{Scod}$ and${I}_{L}$ are functions which take into account the excess of substrates, and the self-shading phenomenon, respectively. They are expressed as:

where ${K}_{Scod}$ and ${K}_{L}$ are the COD and the light inhibition constants.

${I}_{P}$ is the PHB consumption function occurring under substrate starvation, which accounts for the PHB concentration decrease in PNSB cells occurring when organic substrates are in low concentrations. It reads as:

${I}_{N}$is the nitrogen inhibition function, allowing for the switch between the exponential growth phase and the stationary phase. Differently from the other inhibition functions, a secondary substrate inhibition function has been used as reported in Eq (15):

where ${K}_{N}$is the nitrogen inhibition constant.

Simulations were carried out using Matlab (MATLAB R2020b, The MathWorks Inc., Natick, MA). The ODEs system was solved by ODE15s tool. The model has been verified to allow for the COD balance closure. Therefore, the amounts of substrates, biomass, hydrogen, PHB, and other catabolites have been properly reported in mgCOD L^-1. Nitrogen unit has been expressed as mgN L^-1.

All the model outputs have been finally converted in standard units consistently with experimental data. Unit conversions have been calculated for substrate, biomass, hydrogen, and PHB, using theoretical COD values (mgCOD mg^-1), estimated using Eq (16):

The model output of biomass represents the active biomass. Active biomass has been converted to total suspended solid (TSS) content, gTSS L^-1, using the Rhodopseudomonas palustris chemical composition CH_1.8O_0.38N_0.18 reported in ^[34]. The total biomass has been calculated accounting for the PHB concentration incorporated in the microbial biomass. The PHB theoretical COD value has been calculated using its monomeric formula C₄H₆O₂. The PHB output has been successively expressed as a percentage of the total biomass. Hydrogen production has been expressed in terms of volume (mL L^-1) under 1 atm and 25 ℃.

5.   Calibration and sensitivity analysis

The model has been calibrated using experimental data published by Chen et al. ^[28]. The experiment considered for the model calibration has been conducted using acetate as carbon source and glutamate as nitrogen source. Lab-scale serum bottles (200 mL total volume) were fed with 100 mL of culture medium and inoculated with Rhodopseudomonas palustris WP3-5. Reactors were flushed with argon gas and then incubated at 32 ℃, under continuous illumination (6000 lux). Biomass growth, acetic acid concentration, hydrogen production, PHB content, and soluble COD were periodically analysed.

Based on the selected data set, the initial conditions of the modelled variables were fixed at:

${X}_{0} = 261$ mgCOD L^-1;

${S}_{\mathrm{1, 0}} = 2703.79$ mgCOD L^-1;

${S}_{\mathrm{2, 0}} = 0$;

${S}_{\mathrm{3, 0}} = 0$;

${N}_{0} = 4.672$ mgN L^-1;

${PHB}_{0} = 3.74$ mgCOD L^-1;

${H}_{\mathrm{2, 0}} = 0$;

${Cat}_{0} = 0$.

Concerning the model parameters, Table 1 reports all the adopted coefficients in this study. N_bac was calculated as the percentage of nitrogen over biomass from the experimental data used for the calibration. Based on PF modeling literature, the parameters whose value varied in a broad range or were completely absent in other studies were calibrated using a trial-and-error method to fit experimental data. In detail, the model was run several times by modifying each value, until the model properly reproduced experimental data. The selected parameters were calibrated individually. Once the first parameter was calibrated, another parameter was allowed to vary. The algorithm ended when a reasonable fitness with experimental data was reached ^[35]. This method is highly recommended for the calibration of mathematical models with parameters varying in a wide range ^[36]. The order of parameters for the calibration procedure was selected based on a Local Sensitivity Analysis (LSA) study. In detail, parameters with higher sensitivity indexes where calibrated first. LSA is able to investigate the effect of a single parameter variation on a specific model output ^[37]. The influence of a varying parameter is studied by computing the output variable (named quantity of interest) first-order partial derivative over time with respect to the considered parameter. For each parameter, the maximum or the minimum value of the obtained partial derivatives was used as Sensitivity Index (SI), depending on which one is higher in absolute terms ^[38]. The sensitivity Index was estimated as follows:

where ${SI}_{i, j}$ is the ${j}^{th}$ parameter sensitivity with respect to the ${i}^{th}$ output state variable at time $t$ and ${W}_{i, j}$ is the value of the partial derivative of the ${i}^{th}$ output state variable with respect to the ${j}^{th}$ parameter at the instant time $t$. ${C}_{i}$ is the ${i}^{th}$ chosen output state variable concentration, $C$ is the vector of the state variables concentration, ${u}_{j}$ is the ${j}^{th}$ investigated parameter value, $u$ is the vector of the model parameters and $\tau$ is the final instant time.

Hydrogen production and PHB accumulation were set as model output for Eq (17), and the analyzed parameters were: $Y$, ${\nu }_{maxp}$, ${K}_{S}$, ${K}_{P1}$, ${K}_{N}$, ${F}_{1}$, ${F}_{p}$, ${F}_{{cat}_{1}}$ and ${K}_{P}$. The sensitivity analysis was performed by setting the previously mentioned initial condition. The initial values of the investigated parameters were set by varying them arbitrarily to fit experimental data. The Matlab function sens_sys ^[39] coupled with the ODE15s solver were used to perform the numerical analysis. To verify the accuracy of the model calibration, the following performance indicators were determined: the mean absolute error (MAE), its normalized form (NMAE), the modelling efficiency (ME), the index of agreement (IoA) and the fractional mean bias (FB). Such indexes are reported as:

where ${P}_{i}$ and ${O}_{i}$ denote the ${i}^{th}$ predicted value and the ${i}^{th}$ observed value, respectively. N is the number of observed values, $\overline{P}$ and $\overline{O}$ indicates the predicted and observed mean values ^[13].

The aim of the PF process is to maximize the production of hydrogen from a given substrate with a stable biological process. Moreover, PHB accumulation represent a further advantage, allowing for the integrated production of bioenergy and biopolymers ^[40]. Therefore, hydrogen production and PHB accumulation were chosen as model outputs to perform the sensitivity analysis. The results of the sensitivity indexes for each investigated parameter are shown in Figure 1 in log scale.

Results showed that the most sensitive parameters for both hydrogen and PHB outputs were the hydrogen yield on substrate (${F}_{1}$), the catabolites yield on substrate (${F}_{{cat}_{1}}$) and the biomass yield ($Y$). Such results are in accordance with the sensitivity analysis performed by Gadhamshetty et al. ^[19], who developed a kinetic model describing cell growth, substrate consumption, and hydrogen evolution using Luedeking-Piret and logistic equations. In the mentioned work, the most sensitive parameters were found to be the yield coefficient of product formation and the maximum cell concentration. Such analogies remark the crucial role of biomass growth and product yields on the products formation prediction. The main difference between hydrogen and PHB outputs was the role of the hydrogen yield on PHB (${F}_{p}$). As expected, hydrogen production was sensitive to this parameter. Conversely, as PHB is not produced from the consumed PHB, ${F}_{p}$ does not affect PHB accumulation.

The SI values for all model outputs ($X$, $PHB$, ${H}_{2}$, and $Cat$) of kinetic parameters ($Y$, ${\nu }_{{P}_{max}}$, ${K}_{{S}_{1}}$, ${K}_{N}$, and ${K}_{{P}_{1}}$) and model parameters (${F}_{1}$, ${F}_{p}$, ${F}_{{cat}_{1}}$, and ${K}_{P}$) are reported in Figure 2. Regarding kinetic parameters, all model outputs were more sensitive to $Y$ compared to other parameters. However, the most relevant effect of $Y$ was observed on the $X$ variable. ${K}_{N}$ was found to affect hydrogen and catabolites production, as well. However, the effect of ${K}_{N}$ on such outputs is less marked compared to the effect of $Y$. Concerning SI of model parameters (Figure 2(b)), hydrogen and catabolites production were sensitive to ${F}_{1}$ and ${F}_{{cat}_{1}}$, respectively. In particular, ${F}_{1}$ and ${F}_{{cat}_{1}}$ SI were about three times higher compared to $Y$ and, therefore, the most sensitive. Aside from the previously mentioned study by Gadhamshetty et al. ^[19], the present work represents the first application of this procedure to a kinetic model describing the PF process. Therefore, a further comparison with literature is not possible.

Table 1 reports the values of all model parameters. Where possible, a literature value or range has been reported.

Due to the lack of similar models in the literature, only $Y$ and ${K}_{{S}_{1}}$ can be compared with literature ranges. The ${K}_{{S}_{1}}$ value of 18,000 mgCOD L^-1 is in agreement with the literature range reported by Obeid et al. ^[17]. The $Y$ value of 0.85 is the maximum value of the literature range, varying between 0.37 and 0.85 (Table 1). Nevertheless, in the present model, $Y$ was multiplied by the ${I}_{N}$ function, which lowered the biomass yield after nitrogen depletion. From experimental evidence, the first growth phase of Rhodopseudomonas palustris is characterized by a prevalence of anabolic reactions. After nitrogen depletion, biomass does not increase anymore. In this stationary phase, non-growing Rhodopseudomonas palustris increases the products yields from the substrate by shifting from the glyoxylate shunt to the tricarboxylic acid cycle ^[43]. Consequently, the biomass growth represents a crucial aspect in modelling the PF conducted by Rhodopseudomonas palustris. Indeed, the accumulation of microorganisms during different growth phases should be linked to the production rates, related to the specific products of the biological process. Most of the previous reported models described the PF process by modelling the biomass growth, the substrate consumption and the hydrogen production as independent mechanisms ^[16]. To date, the only attempt to define a mathematical model able to link the growth of Rhodopseudomonas palustris and the related product spectrum is represented by Zhang et al. ^[27]. To account for diverse product rates related to different growth phases, these authors developed a mathematical model based on two separate modules. Each module described different growth phases (exponential growth phase and stationary phase) and they were combined using switch functions.

In the present model, the introduction of the ${I}_{N}$ function successfully allowed to describe the switch between the biomass growth phase and the stationary phase, linked to the product formation in both phases. Figure 3 and Table 2 report simulation results of the calibration and performance indicators, respectively.

From Figure 2, it is clear that the model fits properly the experimental data in the different growth phases. In addition, Table 2 indicates that the NMAE and FB error indexes were lower than 5% for $S$ and ${H}_{2}$. The same variables have ME and IoA value greater than 98%. This confirm that the model could efficiently fit the substrate consumption and hydrogen production. Concerning the variables $X$ and $Cat$, ME and IoA indicators were high (>85%). However, the NMAE and FB values were lower than 7%. The negative FB value for the variable $X$ indicates that the model underestimates biomass production. However, from Figure 3 it is possible to notice that the biomass concentration value at the end of the fermentation process was well predicted. Indicators also show that the model does not fit very well nitrogen consumption ($N$). Low performance indicators are due to the great difference between the observed and predicted values at day 2. The result could be due to the indirect nitrogen measure of the observed value, which has been conducted analysing the glutamate concentration using a glutamate assay kit. Probably, the very high glutamate consumption in the medium during the first two days of fermentation did not correspond to the nitrogen assimilation by biomass. Indeed, experimental evidences report that Rhodopseudomonas palustris are able to store nitrogen from glutamate as an intracellular nitrogen source. Consequently, once the glutamate is depleted in the medium, bacteria continue to grow for a short period, using their nitrogen reserve ^[27].

From Figure 3(d), during the first phase of the process, the model underestimated the PHB production, as confirmed by performance indicators. Moreover, the predicted PHB consumption, usually detected at the late stage of photo fermentation process ^[3], was not evident from experimental data. It is worth mentioning that the PHB analytical determination is more difficult and less accurate compared to other observation in laboratory tests. Moreover, the PHB accumulation value are very low and were reported as a percentage of the dry weight. Consequently, the accuracy level of measures and the model fitting further decreased. Of course, this result underline that further research efforts on the mechanisms leading to PHB accumulation in PNSB are still required. Nevertheless, the maximum PHB accumulation has been well predicted by the presented model (Figure 3).

6.   Validation of the model

The model was validated using hydrogen and PHB experimental data reported by Toulopakis et al. ^[29]. From the mentioned study, a data set obtained using acetic acid as substrate and Rhodopseudomonas sp. as microorganism was selected. Based on experimental conditions, the initial conditions of the model were modified as follows:

${X}_{0} = 567$ mgCOD L^-1;

${S}_{\mathrm{1, 0}} = 4248$ mgCOD L^-1;

${S}_{\mathrm{2, 0}} = 0$;

${S}_{\mathrm{3, 0}} = 0$;

${N}_{0} = 9.5$ mgN L^-1;

${PHB}_{0} = 25$ mgCOD L^-1;

${H}_{\mathrm{2, 0}} = 0$;

${Cat}_{0} = 0$.

Moreover, the total duration of experiments decreased from 12 to 8 days.

To verify the accuracy of the model validation, the previously described performance indicators were determined according to Eqs (19)–(23). Specifically, the mean absolute error (MAE), its normalized form (NMAE), the modelling efficiency (ME), the index of agreement (IoA) and the fractional mean bias (FB) were determined for validation purposes.

Figure 4 and Table 3 report simulation results of the validation procedure and validation performance indicators, respectively.

The results showed that the model can predict, with an acceptable fitting, the hydrogen production trend. However, the model slightly overestimates the cumulative production, as the predicted and observed values at the end of the fermentation process were 2737 and 2300 mL L^-1, respectively. The model PHB trend was not in agreement with experimental data. As reported for the calibration results, PHB analytical determination is the most critical in terms of accuracy and PHB production is reported as percentage of the dry weight. However, even in this case, the model predicted very well the maximum PHB accumulation. At the end of the fermentation process, the predicted and observed PHB accumulations were 4.5 and 4.8%, respectively. All considered, such results confirm that the model can be used to predict hydrogen and PHB production.

7.   Additional numerical simulations

As aforementioned, the PF process can be affected by inhibitions mechanisms. Among them, the substrates inhibition and the self-shading phenomenon represent the most significative. Therefore, a numerical study was conducted to test the effect of such mechanisms on the process. In detail, two numerical simulation sets were performed, with different initial concentrations of substrate (S = 2.488, S1/2 = 1.244, S2 = 4.976, S5 = 12.440, S10 = 24.880, S20 = 49.760, and S50 = 124.440 gCOD L^-1) and different initial concentrations of nitrogen (N = 4.76, N1/2 = 2.38, N2 = 9.52, N5 = 23.8, N10 = 47.6, N20 = 95.2, and N50 = 238 mg L^-1). Figure 5 show the results of the simulation conducted with different substrates concentrations.

S represent the substrate concentration of the experimental data used for the calibration. As it can be observed from results, a lower concentration (S1/2) determines the lowering of all products (hydrogen, catabolites and PHB). Slower biomass growth and nitrogen consumption can be observed as well. Conversely, the increase of the substrate concentration (S2 and S5) leads to higher values of all products. In particular, under the concentration S5, the substrate is not completely consumed at the end of the fermentation time of 12 days. Consequently, the consumption of PHB is not observed anymore. As previously remarked, this phenomenon only takes when organic substrates are present at low concentrations. The further increase of the substrate concentration (from S10 to S50) determines a lowering of hydrogen, PHB, and catabolites production, due to the substrate inhibition phenomenon. However, in terms of both hydrogen and PHB generation, the simulation with S10 concentration is still preferred compared to the initial value of S. Due to the high substrate concentration, the S20 scenario leads to a higher PHB concentration with respect to the S condition. On the other hand, hydrogen and catabolites productions are lower. For S20 and S50, the substrate inhibition is evident even in terms of biomass generation and nitrogen consumption. According to numerical simulations, the substrate inhibition takes place when acetate concentration is higher than 12.440 gCOD L^-1. Such value represents the best condition to maximize both hydrogen and PHB.

In Figure 6 the second set of simulation is reported. In this case, simulations were performed by varying the initial nitrogen concentration in bioreactors. The condition indicated with N1/2 represents a nitrogen starvation condition. Indeed, due to the lower biomass concentration, the substrate uptake, the hydrogen, and catabolites productions are quite low. Conversely, PHB accumulation increases. The observed PHB accumulation is due to the lower substrate consumption, leading to a limited PHB depletion phenomenon observed. From N2 to N50, it is possible to observe a progressive inhibition of the whole process, due to the high biomass concentration and the occurrence of the self-shading phenomenon. Moreover, it is possible to observe that the maximum biomass concentration does not exceed the value of about 1200 mgTSS L^-1. Such a value is comparable to results observed in experimental studies ^[29,44]. According to the performed simulations, the best initial nitrogen concentration for hydrogen production is 4.76 mgN L^-1. On the other hand, to increase the PHB accumulation, the lower value of 2.38 mgN L^-1 is advisable. The self-shading phenomenon strongly affects the process when the initial nitrogen concentration is 4.76 mgN L^-1, which correspond to the final biomass concentration of about 600 mgTSS L^-1.

8.   Fast-slow analysis for hydrogen production

The contextual PHB accumulation by PNSB opens an additional question related to the effective maximization of H₂ production from the organic compound ${S}_{1}$. Indeed, the subsequent bioconversion of the accumulated PHB in H₂ leads to a multi-timescale problem where the fast-variables quickly tend towards their steady-states which leads to the stationary phase of the accumulated biomass $X$ and the almost complete nutrient degradation for anabolic reactions. The dynamics of the remaining variables (e.g., the slow-variables) occur on a longer time-scale ^[45]. This is represented by the additional H₂ produced by PNSB from PHB indicating the multi-timescale behaviour. The standard methodology to reduce this system is to make steady-state assumptions on the fast variables and explore the dynamics in the long-time scale ^[46]. To highlight this biological behaviour in the current system, the complete system of Eqs (4)–(9) was considered.

This was first simplified assuming that catabolites $Cat$ were not formed in ideal conditions, and bacteria can exclusively produce PHB and H₂ from the organic substrate. As a consequence, ${F}_{{cat}_{1}} = 0$, ${F}_{p} = 1$, and Eq (9) leads to $\frac{dCat}{dt} = 0$. Based on experimental evidence and two-timescale assumption, the biomass production was supposed to reach a steady condition $X$ = $\widehat{X}$ and Eq (4) was simplified in $\frac{dX}{dt} = 0$, with consequent $\frac{dN}{dt} = -{K}_{d}\widehat{X}$, given by Eq (5). The system reduced to

If nitrogen N→0, as it is usually added in its stoichiometric amount for biomass accumulation (anabolic reactions), the system further reduces to

This system captures the dynamics of the full system on the long-time scale and can be treated as a dynamical system that depends on additional parameters such as the steady-state concentration, $\widehat{X}$ and other determining constants from the fast variables. By adding the latter equations, it is evident that $\frac{d{S}_{1}}{dt}+\frac{dPHB}{dt}+\frac{d{H}_{2}}{dt} = 0$, where ${S}_{1}$, $PHB$, and ${H}_{2}$ are general positive functions for biological meanings, as they represent the concentrations of these components. The system is further reduced in two equations by using ${S}_{1}$ = $-\left(PHB+{H}_{2}\right)+c$, where $c$ is a positive constant $c > PHB+{H}_{2}$. The reduced system was further studied numerically by using the fermentation results of the complete model after 9 d (Figure 3) as initial condition for biomass and PHB concentration. Figure 7 shows simulation results in terms of H₂ production and PHB degradation over time.

The benefit of the fast-slow analysis is the tractability of the reduced system. It is well known that dynamical systems are substantially more well-behaved in two dimensions ^[46], hence the numerical simulations are in some sense generic representations of the dynamics of the full system. From numerical results, the complete PHB is reached after 10 days of simulation, which is technically not feasible as the estimated increase of hydrogen production is around 1%. Moreover, the accumulated PHB can be further extracted to produce bioplastics from microbial cells, leading to an additional value-added route for the exhausted photo fermentative biomass.

9.   Conclusions

A novel mathematical model has been calibrated and validated on experimental data reporting hydrogen and PHB production by Rhodopseudomonas palustis. Results demonstrated that the proposed system of ODEs, accounting for the main process mechanisms and limiting/inhibition factors, is suitable to effectively describe the process. Among model parameters, the yields represent the most sensitive ones and have to be carefully calibrated. From numerical simulations, it emerged that the model is consistent to simulate the substrate and self-shading inhibition phenomena. Other environmental effects can be simply added for future applications, by including other inhibition functions to the current structure of the model. The present model describes the process from a biochemical point of view, but it can be modified to reproduce continuous reactors by adding information related to influent/effluent volumes and components concentrations for designing purposes. It allows for the prediction of the maximum hydrogen and/or PHB production without the need of costly and time-consuming experiments. Further perspectives include the modification of the model to extend it to more complex case studies, such as the PF of waste substrates by mixed phototrophic cultures.

Acknowledgments

Authors would like to acknowledge the project "METAGRO-bioMETanazione dei sottoprodotti della filiera AGROindustriale campana" CUP: B18H19005240009. N. Cogan acknowledges funding from the National Science Foundation (CBET 2210992). This work has been performed under the auspices of the GNFM of INdAM.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.