An Identifiable State Model To Describe Light Intensity Influence on Microalgae Growth.

Despite the high potential as feedstock for the production of fuels and chemicals, the industrial cultivation of microalgae still exhibits many issues. Yield in microalgae cultivation systems is limited by the solar energy that can be harvested. The availability of reliable models representing key phenomena affecting algae growth may help designing and optimizing effective production systems at an industrial level. In this work the complex influence of different light regimes on seawater alga Nannochloropsis salina growth is represented by first principles models. Experimental data such as in vivo fluorescence measurements are employed to develop the model. The proposed model allows description of all growth curves and fluorescence data in a reliable way. The model structure is assessed and modified in order to guarantee the model identifiability and the estimation of its parametric set in a robust and reliable way.

Introduction

Microalgae-based processes are considered one of the most promising alternative technologies for the production of liquid fuels in the transport sector.[1] The main advantages of microalgae with respect to other possible feedstock are the high potential productivity and the absence of competition with traditional crops for arable land and clean water. However, this potential is still theoretical, and algae production on large scale is not profitable yet. Several issues need to be addressed to reach this objective, ranging from algae cultivation and harvesting as well as products extraction.[1,2]

The availability of reliable models would be greatly beneficial as the possibility to represent the fundamental physical, chemical, and biological phenomena would allow one to assess the interactions between equipment design and product yields and to scale-up and optimize the process design and operation.[3]

However, to be reliable a predictive model requires a specific set of parameters to be estimated in the most precise and accurate way. One important advantage of a model whose parameters have been reliably identified is that it can be used to predict the system response also in conditions (significantly) different from those tested during the identification experiments, which is the typical case when a model is exploited for process scale-up and optimization. To achieve that, however, model parameters must be identifiable, and this may present critical issues even in relatively simple models.[4]

In this work we will discuss a modeling structure to represent some key phenomena in algae growth by giving special attention to the identifiability of the proposed model. Algae growth is affected by several variables such as nutrient availability, temperature, and mixing, etc. However, being algae photosynthetic organisms, light is the key variable determining growth efficiency and kinetics: for this reason, the focus of this work will be on the representation of its influence on growth. It is worth stressing that the correlation between light intensity and growth is a very complex one, and while low irradiation is limiting, its excess drives to the formation of reactive oxygen species and has an inhibitory effect.[5]

In the literature several models of photosynthetic biomass growth have been presented. It is possible to divide such models into two main groups: “physiological models” and “state models”. Physiological models attempt to describe the dynamic behavior of photosynthetic cells and propose approximations for the actual mechanisms involved in the cells’ growth. These models may try to represent the optimal allocation of energy and nutrients during cells’ activities (e.g., Ross and Geider[6]) or to represent a specific metabolic reaction (e.g., Marshall et al.,[7] where the damage and repair cycle of protein D1 is described). Usually these models are extremely detailed and involve a large amount of variables and parameters. The actual identification procedure may be extremely complex (sometimes even impossible) and require numerous, highly specific, and costly experiments.

State models are instead based on the concept of a photosynthetic unit (PSU) and are more instrumental for simulating and optimizing industrial cultivation systems. The PSU is defined as the sum of the light harvesting complex, the reaction center, and the associated apparatus, which are activated by a given amount of light energy to produce a certain amount of photoproduct.[8] These models are called state models because the PSU can be in different states of excitation. One of the first proposed state models was the model by Fasham and Platt,[9] whose objective was to describe photoinhibition. Several other modeling approaches have been proposed over the years.[10−16]

In this work we will consider the models by Camacho Rubio et al.[10] and by Eilers and Peeters (initially developed by Eilers and Peeters[11] and then improved by Wu and Merchuk[16]). These models are capable of representing the key phenomena of interest in this work, and they are reasonably simple so as to limit possible identifiability issues. In the Eilers and Peeters model, the authors assume that if an activated PSU absorbs an additional photon, it may become inhibited. For this reason they assume the rate of photoihibition to be proportional to light intensity. It is also assumed that photosynthesis (and by consequence biomass growth) is proportional to the transition between the activated state and the resting state. Later Wu and Merchuck[16] modify the model of Eilers and Peeters, introducing a constant maintenance factor in the description of biomass growth.

In the model by Camacho Rubio et al.[10] both photoinhibition and photoacclimation, which represent the cells’ response to optimize their photosynthetic apparatus to different light intensities, are considered (note that in the original work, and in several others, photoacclimation is instead called photoadaptation: photoacclimation is however a more accurate definition: in fact, adaptation refers to the organisms’ modification during evolution to their environment, thus responses with time scales extremely longer that the ones considered here). Photoacclimation was at first represented as a steady state process but more recently extended by the same authors in order to represent its dynamics, together with the effect of nonphotochemical quenching (NPQ) and dark respiration (Camacho Rubio et al.[12]). This last model is indeed very flexible, but the number of model parameters to be estimated is very high and their precise identification may become a long and difficult task, especially if a limited amount of data is available, and for this reason it will not be discussed further in this work.

For model development and identification we consider experimental data referring to a particular species of microalgae of industrial interest (Nannochloropsis salina (N. salina)) grown in nonlimiting nutrients conditions and in a flat-plate photobioreactor.[17] These data sets were selected because experimental conditions were optimized to minimize all influences on algae growth other than light intensity. In fact, nutrients and CO2 were provided in excess, but also the photobioreactor light path was minimized to reduce as much as possible light attenuation due to cells’ shading and scattering. Sforza et al.[17] demonstrated that this assumption was an acceptable approximation, especially considering that we are interested in representing the exponential growth phase, where nutrient availability is high and cell concentration low. Accordingly, these data represent an accurate description of the influence of light alone on algae growth, minimizing the effect of other parameters. It is worth clarifying that other phenomena, which also play a major influence on algae productivity in industrial photobioreactors, such as the dark/light cycles due to mixing are not considered in this model.

There are two main contributions in the present work. On the methodological side, a step by step approach will be presented and applied to guarantee the identifiability of the growth model. Second, the utilization of both biomass concentration (growth curves) and multiple fluorescence measurements will allow shading light on some fundamental phenomena in the correlation between illumination and growth; in particular, differently from previous contributions, measurements of the light profile of PSU saturation will be exploited.

The work is organized as follows: after outlining the general identification methodology, the available experimental data and the two candidate models will be introduced and discussed. The successive section is about model discrimination and the enhancement of the selected model. Then, an identifiability analysis and a reparameterization approach will allow setting up an identifiable model. The performance of the model in describing algal growth will be critically discussed. Some final remarks will conclude the work.

Model Developing Approach

Identifiability is a key issue to guarantee reliability and predictive capability in a model being developed. Figure 1 outlines the basic tasks and information flux required to achieve such a target. The preliminary step is to identify the phenomena that need describing and the fundamental mathematical laws that should be implemented to represent them. Here we assume that some modeling assumptions are already available. In other cases, preliminary experimental data may be needed to envisage the correlation among data and to set up a suitable physical interpretation through a mathematical model.

Figure 1

Information flux of the model identification procedure.

Typically, some available data may be exploited at this stage to choose among competitive modeling approaches through suitable discrimination techniques.[18,19] At least in the easiest cases a χ2 test on experimental data may be sufficient to make the discrimination.[20,21,41]Ad hoc experiments can also specifically be designed to allow for a more effective and reliable discrimination among different candidates.[21−22] Once a suitable candidate model has been selected (and a preliminary estimation of its parameter has been carried out), the model may need upgrading to improve its capability of representing the phenomena being investigated (new experiments may be needed and possibly designed, and an estimation of all model parameters should be attempted).

Then it is extremely important to verify the model identifiability, i.e., to confirm that the optimal set of parameters values is unique and that their values can be determined in a precise way (and ideally in a physically meaningful way). If some identifiability issues arise, then a first approach to tackle the problem is to reparameterize the model.[23] If this is not enough, then the model structure should be modified.

Once the model is proved to be identifiable, the final parameter estimation can be performed. Note that even when a model is identifiable, measurements of noise and other uncertainty effects may still hinder its practical identifiability, although properly designed experiments may help in tackling the issue.[24] Once the final parameters estimation has been performed, new data, not involved in model calibration, should be used to validate the model. This approach has been applied to the specific case study and is discussed in the following sections. For simulation and parameter estimation purposes, the gPROMS software has been used.[25]

Modeling Approaches

The main advantage of state models is that they reduce the complexity of photosynthesis into a few possible states of the PSUs. This simple structure is also particularly effective in the use of fluorescence measurements, which can be exploited to monitor the PSU populations at different states. A PSU is defined as the sum of the light harvesting complex, the reaction center, and the associated apparatus, which are activated by a given amount of light energy to produce certain quantities of photoproducts. The Eilers and Peeters model in the form proposed by Wu and Merchuck[16] (afterward denoted as EPM) and the Camacho Rubio model (Camacho Rubio et al.,[10] later called CRM) are two of the simplest models that can describe photosynthetic biomass growth as a function of light intensity. Both models consider that a PSU can assume three different states of excitation: (1) the resting (or open) state, which is the state of PSU before the light energy excites the reaction center; (2) the activated (or closed) state, that is the state of PSU excited by light energy; (3) the inhibited state, which is the state of PSU damaged by an excess of light energy. Both models do not consider any limitation on nutrients availability or mass transport of nutrients: i.e., only the exponential growth phase is described. CRM considers photoacclimation too, but without any representation of its dynamics. This is a reasonable assumption also in our case study, since acclimation characteristic time scale (hours or days) is significantly larger than the time scales of the other phenomena being investigated.

In the work of Wu and Merchuck[16] EPM was used to fit the data of experiments carried out in a thin tubular loop reactor. Part of the reactor was kept in dark to simulate mixing; light intensities used for the experiments were 110, 220, and 550 μE/(m2 s). The measurements used to calibrate the model were both biomass concentration and dark fluorescence measurements (fluorescence measurements will be described in section 4). Each experiment was carried out for 48 h, and measurements were taken every 12 h. In the work of Camacho Rubio et al.,[10] CRM was applied to a wider range of light intensities (ranging from 0 to 2000 μE/(m2 s)) and to different light regimes (constant light, flashing light, and day–night cycle). Data used by the authors were growth rate constant and P–I curves taken from the literature.

EPM assumes that the number of PSUs is constant with respect to light intensity and accordingly refers to the PSU x1, x2, and x3 to represent the resting, activated, and inhibited states, respectively. Conversely, CRM assumes that the number of PSUs is a function of light intensity (indicated as at), and the model equations are expressed as a function of the amount of the PSUs in each of the three states (a1, a2, and a3). Parts a and b of Figure 2 illustrate the two model structures.

Figure 2

(a) Scheme of EPM. x1, x2, and x3 represent the fractions of PSUs in resting, activated, and inhibited states, respectively. (b) Scheme of CRM. a1, a2, and a3 are the numbers of PSUs in resting, activated, and inhibited states, respectively; at represents the total number of PSUs.

Eilers–Peeters model

In Figure 2a, the three PSU states are represented by circles and the possible state transitions are represented by the arrows. The resting state PSU can capture light energy and transfer it to an activated state. The PSU in the activated state can be damaged by light, or pass down the energy to start the dark phase of photosynthesis (and then return to a resting state). An inhibited PSU can be recovered and then return to the resting state. The reaction rate of the transitions involving the absorption of light (i.e., x1 → x2 and x2 → x3) is assumed to be first order with respect to light intensity. The other two transitions are assumed to be zero order with respect to light intensity. Each transition is assumed to be first order with respect to the PSU fraction involved in the transition. The growth rate constant (μEP (h–1)) is assumed to be proportional to the state transition from activated to resting state, representing the photochemical reactions. Considering that the growth rate can be negative in the dark or at very low light intensity, a constant maintenance (MEP (h–1)) factor is introduced. The model equations are as follows:The set of parameters (whose physical meaning is summarized in Table 1) is represented by vector θ̂EP = [kaEP,kdEP,kiEP,krEP,kpEP,MEP].

Table 1

Parameters of EPM Significance and Units

param	significance	units
kaEP	kinetic constant of the activation reaction rate	m2/μE
kdEP	kinetic constant of the deactivation reaction rate (photochemical quenching)	s–1
kiEP	kinetic constant of inhibition reaction rate	m2/μE
krEP	kinetic constant of the recovery reaction rate	s–1
kpEP	proportionality factor between photochemical quenching and biomass growth rate constant	s/h
MEP	maintenance factor	h–1

Camacho Rubio Model

In CRM, the photoinhibition rate is assumed to be proportional to the sum of resting state and activated PSUs, i.e., the active PSUs represented by a1 + a2 [PSUs/cells] in Figure 2b. As in EPM, PSU activation reaction is assumed to be first order with respect to light intensity and to the amount of resting state PSUs. However, in CRM the transition from activated to resting state, related to biomass growth, is defined as a Michaelis–Menten kinetic, assuming an enzymatic reaction as the limiting step of this process. Also, from an analysis of experimental data, the authors[10] assume that the photoinhibition reaction rate is first order with respect to the square root of light intensity. As in EPM the recovery of damaged PSUs is assumed to be a first order reaction with respect to the number of damaged PSUs. Finally, as anticipated, CRM includes photoacclimation. The total amount of PSUs (at [PSUs/cells]) in CRM is thus assumed to be a hyperbolic decreasing function of light intensity. As for EPM, the growth rate constant (μCR (h–1)) is assumed to be proportional to the transition from activated to resting state and a constant maintenance factor (MCR (h–1)) is introduced. The model equations are as follows:As reported in the literature, it appears that the condition KSCR ≫ a2 is true for the entire range of light intensities considered in the experiments (ranging from 50 to 1000 μE/(m2 s)). Thus, the Michaelis–Menten kinetic may be well approximated by a first order kinetic as in EPM. Accordingly, the reduced set of parameters is represented by vector θ̂CR = {kaCR,kdCR,kiCR,krCR,kpCR,kcCR,MCR} (parameters’ physical meanings and units are reported in Table 2).

Table 2

Parameters of CRM Significance and Units

param	significance	units
kaCR	kinetic constant of the activation reaction rate	m2/μE
kdCR	kinetic constant of the deactivation reaction rate (photochemical quenching)	s–1
kiCR	kinetic constant of inhibition reaction rate	m/(μE·s)0.5
krCR	kinetic constant of the recovery reaction rate	s–1
kpCR	proportionality factor between photochemical quenching and biomass growth rate constant	s/h
kcCR	rate constant involved in the photoacclimation process
MCR	maintenance factor	h–1

Experimental Setup and Available Data

The aim of this work is to describe the growth of microalgae in nonlimiting nutrient conditions and according to the hypothesis that the light intensity is constant with respect to the culture time and depth. The fundamental phenomena to be described are as follows: (i) reaction centers oxidation/reduction cycle, to represent the photosynthesis, and (ii) the damaging effect of excess light on PSUs (photoinhibition).

Our data refer to algal cultures grown at different light intensities.[17] During the experiments, microalgae were acclimated to the light used and to the geometry of the photobioreactor. Each experiment was conducted in parallel at least twice and in two identical photobioreactors, in order to ensure its reproducibility.[17] Only the data in the exponential phase of the original growth curves were used for the parameters estimation, since our model represents only exponential growth and does not consider nutrients limitation. Experimental setup was built to limit as much as possible cells shading by decreasing as much as possible the light path and working at low cells concentration. This, together with the presence of nutrients and CO2 in nonlimiting amounts, ensures that growth is dependent only from the light intensity reaching the culture.

Moreover, experimental measurements of fluorescence will be used as additional data for the parameters estimation.[26] These data are commonly available for the photosynthetic organisms and have been exploited to estimate photosynthetic efficiency in a large body of experimental literature (reviewed by Maxwell and Johnson[26]). In our case we considered in particular two fluorescence parameters, related to the oxidation state of the algal cells: the parameters q (or Fv/Fm),[27−30] and qL, whose physical meanings will be briefly detailed in the following.

Parameter v/Fm (with Fv = Fm – F0) was measured for each light intensity at which microalgae were grown. F0 represents the fluorescence in the dark when all reaction centers are open. Conversely, the maximum Fm is the fluorescence after a saturating light flash, when all reaction centers are closed. Here the fluorescence is higher than F0 because saturated PSUs are unable to perform photochemistry and therefore a larger amount of excitation energy is re-emitted as fluorescence. When PSUs are active, there is a large difference between Fm and F0 (and thus a larger Fv/Fm), while if PSUs are photoinhibited, Fm is closer to F0. For this reason, the Fv/Fm parameter is commonly used in the literature to quantify active PSUs in vivo. The precise value of Fv/Fm in the case of fully active PSUs is known to be variable between species, since it depends on specific properties such as the antenna size. Here it is set to be equal to 0.65 as this is a value commonly measured in several microalgae[31] and also in healthy Nannochloropsis cultures, exposed to low light.[17,30] In order to have a parameter bound between 0 and 1, it is possible to define the parameter qnorm = q/qmax, where qmax is the value of Fv/Fm in the case of fully active PSUs.

In cultures exposed to different light intensities, a decrease in this value indicates the presence of photoinhibited PSUs, as normally experienced in high light conditions. For this reason, it can be used to estimate the content of photoinhibited PSUs and thus it can be correlated to x3 and a3 populations, using the definition of EPM and CRM, respectively.

While q is rather commonly employed in several similar studies,[16] in order to have a better representation of the oxidation state of the PSUs in illuminated cells, here we also included the fluorescence parameter qL, which provides a linear estimation of the saturation level of PSU as discussed in detail in the work by Kramer et al.[32] This means that qL is 1 when all active PSUs are open, while it decreases to 0 when all PSUs are closed and photosynthesis is saturated. Parameter qL = [(Fm′ – FS′)/(Fm′ – F0′)](F0′/FS′) was measured for 21 different light intensities with a PAM fluorometer. Here Fm′ and F0′ stands for the same minimal and maximal fluorescence, measured in light adapted cells, while FS′ stands for stationary fluorescence in illuminated cells. Thus, qL can be exploited as an estimation of the relative ratio of x1 and x2 populations (a1 and a2 according to CRM). The complete set of experimental measurements used in this work is reported in Appendix A.

Model Discrimination and Preliminary Parameters Estimation

The first objective is to discriminate between the two alternative models so as to select the most suitable one to describe our system. Parameter estimations were performed for both models, based on the entire set of experimental data. The different performance is quantitatively summarized by the χ2 test (in the case of EPM, we have χ2 = 890.1, whereas, for CRM, χ2 = 301.6).

Figure 3 shows the behavior of the two models in representing fluorescence experimental profiles. Although both fittings are quite unsatisfactory, CRM clearly outperforms EPM. In fact, Figure 3a shows that although EPM provides a slightly better fit for the profile of q at low light intensities, only CRM is capable of representing the trend at a higher intensity. Figure 3b, however, shows that although EPM performance is still the worst one, both models cannot properly represent the oxidative state of PSU in light adapted cells, consistent with the fact that these kinds of measurements were not considered in building such models. This means that both models are not accurate in estimating the PSU oxidative state.

Figure 3

Measurement (black circles) and predicted values of (a) qnorm and (b) qL. Red solid lines represent the profiles according to EPM, while the dashed green lines represent the profiles according to CRM.

Also in the predictions of the growth rate constant, reported in Figure 4, CRM outperforms EPM. In fact, EPM does not predict the correct value of light intensity at which the maximum growth rate is reached and underestimates the growth rate constant at low light. CRM correctly predicts the optimal light intensity but it overestimates the growth rate constant at low light.

Figure 4

Growth rate constant predicted by the EP (solid line) and CR models (dashed line) and experimental values of the growth rate constant (black squares).

No additional experiments are needed for discrimination purposes and CRM is then selected, although the experimental data demonstrate that further improvements are required to improve the fitting in the region between 150 and 700 μE/(m2 s) with respect to the PSU oxidation state. It is worth underlining that the main difference between EPM and CRM is the presence, in the latter, of a photoacclimation term expressed through the number of PSUs per cell. This suggests that the inclusion of this response is absolutely necessary for an accurate description of the photosynthetic performances and is consistent with its biological relevance, demonstrated by the fact that acclimation responses are conserved in all photosynthetic organisms.

Enhancing CRM

To improve the model, a more detailed description of some fundamental biological phenomena needs to be introduced. According to several works in the literature, PSII photoinhibition occurs at all light intensities.[28,33,34] Therefore, we assumed that photoinhibition does not depend on the number of active PSUs, but is simply related to light intensity. Above Icr (the light intensity where photosynthesis is saturated) a second process (photoprotection) is activated: in these conditions a fraction of the energy absorbed does not lead either to photochemistry or to PSU damage but is simply dissipated (e.g., because of NPQ). Accordingly, eq 6 has been modified to represent the different behavior above and below Icr. Furthermore, in eq 6, the reaction order with respect to the light intensity is assumed to be 0.5: since there is no clear physical reason for setting such a value, here we decided to increase the model flexibility and its capability of incorporating all energy dissipation phenomena, by making the reaction order a parameter (α) to be estimated.

As mentioned above, photoprotection mechanisms are activated when photochemical reactions are saturated, and accordingly parameter Icr plays a key role in representing this behavior. The value of critical light intensity when photosynthesis is saturated has been fixed to 150 μE/(m2 s), which is the light intensity in response to which Nannochloropsis growth is maximal, and the limit over which the growth rate is not linearly dependent on light anymore.[17] This choice is well consistent with the observation of NPQ dependence from light intensity observed in Nannochloropsis cells, which shows activation only over this limit (NPQ data are reported in Appendix B). This parameter thus depends on both light intensity and the number of active PSUs.

Also eq 8 representing photoacclimation is modified to allow for a higher flexibility: the hyperbolic form is retained, but the light exponent becomes an additional parameter (α2) to be estimated (see later on, eq 13).

Finally, the representation of the maintenance factor is modified, too. In CRM the maintenance factor is treated as a constant (in fact, this is a typical assumption). However, the maintenance factor should vary with light intensity to account for the metabolic cost of repairing damaged PSUs.[35] In order to do this, the easiest way is to express the maintenance factor as a linear function of the damaged PSUs fraction. Since the fluorescence measurements return the number of active PSUs, the variable term of the maintenance factor was related to the difference between the maximum value of fluorescence (qmax) and the current value of fluorescence (q). Finally the fluorescence measurements q and qL are related to the oxidation state of the PSUs, as discussed in section 4.

The modified CRM (mCRM) is thus constituted by the following set of equations:where a new vector of the model parameters is θ̂ = [ka,kd,ki,0,ki,1,kr,kp,kc,kM,M0,α,α2].

Identifiability Analysis

In order to be reliable and suitable for process simulation and optimization, the model for photosynthetic biomass growth has to be identified against the experimental data. It can be verified that mCRM shows identifiability issues if a parameter estimation is performed on the whole parameters vector θ̂. The estimation of the model parameters is characterized by large confidence intervals for some parameters, and more than one set of optimal values can be determined (i.e., the model is not uniquely identifiable). In order to overcome this problem both global (or structural) and practical identifiability of the model have been studied. First of all, the global identifiability has been verified using a differential algebra based method. Afterward, the practical identifiability has been assessed through a sensitivity analysis and a model reparameterization has been performed.

Global Identifiability Analysis

The first step for testing model identifiability is represented by global identifiability analysis. Global identifiability analysis is in fact a necessary condition for practical identifiability and can provide the minimum number of observations required to identify an unique set of optimal parameter values. The two hypotheses, upon which structural identifiability analysis rely, are as follows: (i) complete absence of measurement errors and (ii) a perfectly accurate model structure. Those two assumptions refer to an ideal case, and therefore, once a model is verified to be globally identifiable, practical identifiability, too, needs assessing. Several methods to study the structural identifiability of nonlinear models have been developed in the past two decades and are nicely reviewed in the work by Miao et al.[36]

A rather common approach is given by the series expansion method (e.g., Dochain et al.[37] applied this approach to kinetic models of activated sludge respiration). It requires that functions representing the model are infinitely differentiable, as the method involves the calculation of arbitrary order derivatives. Several examples of application of this method can be found in the literature. However, the series expansion method has a serious drawback: for high dimensional models, high order derivatives are necessary and the resulting equations can easily become too complicated to solve. In fact, it was verified that in this case the series expansion method leads to an intractable problem. Thus, in order to overcome the difficulties related to high order derivatives calculation and the resolution of the resulting equations, a method based on differential algebra was taken into account. In particular, the software package DAISY[38] has been used here. DAISY, as with all differential algebra based methods, requires equations of the models to be written in the form of differential polynomials.[36] For this reason, mCRM has been slightly modified as detailed in Appendix C.

Results indicate that a single experiment is not sufficient to identify a unique value of all of the parameters. However, if (at least) three parallel experiments are carried out at different light intensities and both concentration and fluorescence parameters are measured during the experiments, the model is globally identifiable. This means that the experienced identifiability issues depend on practical identifiability. To tackle the problem, a sensitivity analysis and a reparamereterization will be performed. Note that in several cases, where more complex models need considering, the analysis of the model practical identifiability may be the only test that can be carried out, since the verification of global identifiability may not be viable.

Sensitivity and Correlation Analysis

As the model is structurally identifiable, its practical identifiability will now be tested through a sensitivity analysis. The sensitivity q of the ith response to the kth model parameter is defined aswhere y is the ith measured responses predicted by the model, y′ is the same response obtained from a perturbed value of the kth parameter θ, and Δθ is the perturbation (NM is the number of measured responses and Nθ is the number of model parameters). The principal goal of sensitivity analysis is to evaluate the impact of each parameter on the measured responses and to underline the presence of correlation among specific subsets of model parameters. In an overparameterized model, near-zero sensitivity values would be obtained, leading to the nonidentifiability of some subsets of model parameters. In the most desirable case, sensitivity profiles should be clearly distinct and far from being symmetrical (i.e., they should present a low mutual correlation). As the sensitivity analysis requires prespecified parameter values, values by a preliminary parameter estimation have been considered. Perturbation Δθ was set equal to 1% of the parameter values.

First, the sensitivity profiles of biomass concentration will be taken into account. Although not shown here for the sake of conciseness, it was verified that all sensitivity profiles behave as exponential curves. As a consequence, the final value of sensitivities (value at 120 h) is sufficient to analyze the system behavior in different experimental conditions. The sensitivity profiles evaluated at four different light intensities are reported in Figure 5.

Figure 5

Final values of dynamic sensitivities for mCRM, evaluated at different light intensities.

For light intensities under the critical value, the sensitivities related to parameters ki,0 and kr show opposite values. Above the critical value of light intensity, the sensitivity of kr is the opposite of the sum of the sensitivities of ki,0 and ki,1. This suggests that it may be difficult to exploit biomass concentration measurements to identify ki,0 and kr.

Another critical aspect is related to the fact that parameters ka and kp exhibit a very similar sensitivity at low light conditions (availability of growth data at high light intensities may be necessary for a robust estimation of these parameters). Also it should be noticed that the parameter related to photoacclimation, kc, and the parameter representing the maintenance factor in the dark, M0, exhibit a very similar (and low) sensitivity in all light conditions.

Fluorescence profiles (q and qL) are also considered. Because the dynamics of PSUs are fast with respect to the sampling time, and our measurements are steady state measurements, only the steady state values of sensitivities will be reported. In Figure 6a,b the values of the sensitivities of q and qL are reported for four different light conditions. The sensitivities of parameters α2, kc, kp, kM, and M0 are not reported, since they are zero for both fluorescence measurements (they are related to biomass growth and are not concerned with the oxidation state of the PSUs).

Figure 6

Final values for dynamic sensitivities for the Camacho Rubio model evaluated at different light intensities considering literature values for model parameters.

Considering the steady state sensitivity values of the fluorescence measurements, the following critical aspects can be noticed: (1) parameters ka and kd do not affect the dark fluorescence measurements, while they are completely anticorrelated if light fluorescence measurements are available, showing opposite sensitivities in all light conditions; (2) parameters ki,0, ki,1, kr, and α are not affected by light fluorescence measurements; (3) parameter kr is anticorrelated with ki,0, under the critical value of light intensity, and with ki,1, above the critical value of light intensity; (4) the sensitivity of ki,0 is always quite small.

The results suggest that a model reparameterization may be necessary to help tackling the issue.[23]

Reparameterization of the Camacho Rubio Model

Because the fluorescence data are measurements of the PSU oxidation state under steady state conditions, an analytical expression for fluorescence measurements can be derived from eqs 10 to 12:with R0 = ki,0/kr, R1 = ki,1/kr, and R2 = ka/kd. The interesting aspect is that only four parameters affect the steady state values of the fluorescence measurements: the three ratios R0, R1, and R2 and parameter α. This suggests that, in order to have a practically identifiable model, the values of two parameters affecting the PSU dynamics have to be fixed if no dynamics in the fluorescence data can be included in the data set. It was chosen to fix the values of kr and kd, as those parameters represent the rate constant of the recovery and deexcitation processes, respectively. An approximated estimation can be obtained, considering the time scales of the processes involved: according to literature,[39] values of 100 s–1 for ka and of 2.22 × 10–4 s–1 for kr were assumed. The very low sensitivity to available measurements and the high correlation of kc and M0 make it impossible to identify the two parameters. However, in the case of M0, previous experiments suggest it to be between 5% and 10% of the maximum growth rate.[40] Here we set M0 = 1.5 × 10–3 h–1 (∼8% of the maximum growth rate). For parameter kc the preliminary estimation suggests that a “small” value is required for a good description of data. Since in eq 13 we have that Iα ≫ 1, we verified that kc = 1 is a good approximation for representing all experimental conditions. After reparameterization, the vector of parameters to estimate is θ̂* = [R0,R1,R2,kp,kM,α,α2].

Results and Discussion

In the first part of this section the mCRM parameter estimation results will be presented and discussed. In the second part two additional growth curves and two new measurements of dark fluorescence will be used to validate the model.

Model Calibration

In the estimation procedures, parameters are normalized with respect to the initial values obtained by the preliminary parameter estimation, to increase numerical robustness. The results of parameters estimation are reported in Table 3, along with the confidence intervals and the t-values (for a statistically precise estimation of a model parameter the t-value has to be greater than a reference t-value). The t-value statistic shows that the parameter values are estimated in a statistically satisfactory way.

Table 3

Estimated Values of Parameters of the Reparameterized mCRM, Normalised Values (with Respect to the Initial Values) of the Parameters, Confidence Intervals (conf int), and t-Student Values for Parametersa

param	estd value	normalized value	95% conf int	t-value 95%
R0	4.93 × 10–4	1.12	0.21	5.20
R1	1.34 × 10–2	0.65	0.21	3.07
R2	3.30 × 10–3	2.77	0.41	6.72
kp	1.48 × 10–6	0.41	0.12	3.35
kM	1.12 × 10–1	1.96	0.91	2.15
α	0.45	0.45	0.048	9.52
α2	0.30	0.30	0.078	3.89

The reference t-value is equal to 1.67.

The profiles of fluorescence, predicted after the identification of the reparameterized mCRM, are reported in Figure 7. We can observe that the model correctly fits both the dark fluorescence profile (Figure 7a) and the light fluorescence measurements (Figure 7b). Biomass growth profiles are rather well represented by the model as illustrated in Figure 8, where the six different illuminating conditions are represented.

Figure 7

Measurement (black circles) and predicted values of (a) q and (b) qL. Blue solid lines represent the profiles predicted by the modified Camacho Rubio model. Red stars in panel a represent the experimental data used for model validation.

Figure 8

Biomass concentration profiles at different light intensities predicted by the modified Camacho Rubio model. Black circles represent the experimental measurements.

In Figure 9 the growth rate constant predicted by the model is reported along with the experimental value (i.e., the value obtained fitting the experimental data of growth during the exponential phase with an exponential curve). We can observe that, for all light intensities at which an experiment was carried out, the model well describes the experimental values of the growth rate constant. This suggests that, thanks to the fundamental input of the fluorescence parameters in illuminated cells (qL), the model is capable of reproducing with sufficient accuracy the basic processes of photosynthesis: photochemistry, light damage, and also energy dissipation.

Figure 9

Growth rate constant predicted by the modified Camacho Rubio model (solid line) and experimental values of the growth rate constant (black circles). Red stars represent the experimental values of the growth rate constant for the experiments used in the model validation.

From a statistical point of view, the predicted profiles have a χ2 value of 117.7 (whereas χ2 = 890.1 in the case of EPM and χ2 = 301.6 for CRM).

Model Validation

In order to validate the model, two experiments, not used for model calibration, will be taken into account. In particular two growth curves at 350 and 750 μE/(m2 s) have been considered. For both cultures also the value of q has been measured and was exploited for model validation. In Figures 7a and 9 the validation points are represented by red stars. In Figure 10 the growth curves have been reported.

Figure 10

Biomass concentration profiles at different light intensities predicted by the modified Camacho Rubio model. Red stars represent the experimental measurements used for the model validation.

We can observe that the predictions are accurate for the parameter q. Also the growth curves are predicted in a sufficiently accurate way, although the growth at 350 μE/(m2 s) is slightly overestimated by the model, while the growth at 750 μE/(m2 s) is somewhat underestimated in the final part. From a statistical point of view, validation leads to a χ2 value of 226.23.

Conclusion

A literature model[10] was selected and modified to describe microalgae growth through a rigorous identification procedure. The estimation of the model parameters was performed considering experimental data (growth profiles and fluorescence measurements) on N. salina. Results show that the developed model well represents biomass growth over a wide range of light intensities. The modified model was also able to reproduce fluorescence measurements, including the light profile of PSU saturation. This suggests that the proposed model is accurate enough to represent all major processes of photosynthesis, photochemistry, PSU damage, and energy dissipation. While results in reproducing experimental data are fully satisfactory, it should be underlined that algae growing in an industrial scale photobioreactor are exposed to different conditions. In particular light is not homogenously distributed because of cells shading, and illumination intensity is not constant because of diurnal changes and cells mixing. Finally nutrient availability can also be limiting. Future efforts will be made to expand the model to include these phenomena, also by designing appropriate experiments.

Table 4

Measured Biomass Concentration Profiles at Different Light Intensitiesa

	c (g/L)
time (h)	I = 50 μE/(m2s)	I = 120 μE/(m2s)	I = 150 μE/(m2s)	I = 250 μE/(m2s)	I = 350 μE/(m2s)	I = 550 μE/(m2s)	I = 750 μE/(m2s)	I = 1000 μE/(m2s)
0	9.93 × 10–2	1.81 × 10–1	1.42 × 10–1	5.49 × 10–2	1.30 × 10–1	1.80 × 10–1	1.21 × 10–1	3.02 × 10–1
24	1.44 × 10–1	3.07 × 10–1	2.40 × 10–1	1.10 × 10–1	1.69 × 10–1	–	1.61 × 10–1	4.12 × 10–1
36	–	–	3.60 × 10–1	–		–		–
48	1.57 × 10–1	4.69 × 10–1	4.68 × 10–1	2.21 × 10–1	2.02 × 10–1	4.10 × 10–1	1.86 × 10–1
60	–	–	5.76 × 10–1	–		–		–
72	2.07 × 10–1	8.01 × 10–1	–	3.32 × 10–1	4.06 × 10–1	6.19 × 10–1	3.50 × 10–1	4.56 × 10–1
84	–	–	1.07	–		–		–
96	–	1.43	–	5.57 × 10–1	5.17 × 10–1	7.87 × 10–1	4.80 × 10–1	5.70 × 10–1
120	–	–	1.78	6.92 × 10–1		9.50 × 10–1	6.30 × 10–1	1.22
144	6.31 × 10–1	–	–	8.32 × 10–1		1.61		1.50

Experiments carried out at 350 and 750 μE/(m2s) have been used for mCRM model validation.

Table 5

Values of Growth Rate Constant and Coefficient of Determination R2 Obtained from the Linear Regression of the Growth Curves, Reported in a Semilog Scalea

I [μE/(m2s)]	μ (h-1)	R2
45	0.0126	0.976
120	0.0212	0.998
150	0.0215	0.987
250	0.019	0.954
350	0.0152	0.953
550	0.0144	0.984
750	0.0145	0.975
1000	0.0107	0.875

Data points used for the growth rate constant determination are reported in Table 4.

Table 6

Measured Value of Parameter q and Its Normalized Value qnorm at Different Light Intensity (Normalization Factor Set to 0.65, the Value of q Corresponding to Fully Active PSUs)a

I [μE/(m2s)]	q	qnorm
0	0.650	1.00
46	0.645	0.992
120	0.611	0.940
150	0.592	0.911
250	0.533	0.820
350	0.509	0.783
550	0.503	0.773
750	0.516	0.794
1000	0.499	0.769

The measurements taken at 350 and 750 μE/(m2s) have been used for model validation.

Table 7

Measured Value of Parameter qL at Different Light Intensity

I [μE/(m2s)]	qL
0	0.650
14	0.567
21	0.545
30	0.536
45	0.529
61	0.516
78	0.519
103	0.495
134	0.464
174	0.438
224	0.405
281	0.369
347	0.326
438	0.277
539	0.231
668	0.203
833	0.168
1036	0.138
1295	0.114
1602	0.0871
1960	0.0592