High-Fidelity Modelling Methodology of Light-Limited Photosynthetic Production in Microalgae.

Reliable quantitative description of light-limited growth in microalgae is key to improving the design and operation of industrial production systems. This article shows how the capability to predict photosynthetic processes can benefit from a synergy between mathematical modelling and lab-scale experiments using systematic design of experiment techniques. A model of chlorophyll fluorescence developed by the authors [Nikolaou et al., J Biotechnol 194:91-99, 2015] is used as starting point, whereby the representation of non-photochemical-quenching (NPQ) process is refined for biological consistency. This model spans multiple time scales ranging from milliseconds to hours, thus calling for a combination of various experimental techniques in order to arrive at a sufficiently rich data set and determine statistically meaningful estimates for the model parameters. The methodology is demonstrated for the microalga Nannochloropsis gaditana by combining pulse amplitude modulation (PAM) fluorescence, photosynthesis rate and antenna size measurements. The results show that the calibrated model is capable of accurate quantitative predictions under a wide range of transient light conditions. Moreover, this work provides an experimental validation of the link between fluorescence and photosynthesis-irradiance (PI) curves which had been theoricized.

1 Introduction

Microalgae are considered one of the most promising feedstocks for medium-term biofuel production [1], and have been the subject of extensive research in recent years. Despite this recognised potential, one of the main issues to address is bridging the gap between maximal theoretical biomass productivity (or even lab-scale realised productivity) and the actual biomass productivity in large scale production systems. To meet this objective, models capable of reliable and quantitative prediction of all key phenomena affecting microalgae growth, particularly light utilization, can be a great help to improve our understanding as well as optimise process design and operation [2, 3].

Quite a large number of modelling approaches have been proposed over the years, and it is beyond the scope of this article to provide an exhaustive review of the literature. The focus here is on the so called state models, based on the concept of photosynthetic unit (PSU), which have proved especially useful to characterise photosynthetic operation in practice. Kok [4] first published a state model describing the damage of the photosynthetic apparatus under high light conditions and the subsequent repair processes occurring in dark conditions. Later, Eilers and Peeters [5] proposed a state model that predicts steady-state photosynthetic productivity, also considering the inhibitory effect of excess light. Wu and Merchuck [6] extended Eilers and Peeters’ model by introducing a maintenance term representing the loss of biomass due to respiration when the culture is kept in the dark. The model by Han et al. [7] is equivalent to the model by Eilers and Peeters, but it refines the biological interpretation of the model parameters. In the work by Rubio et al. [8] the effect of photoacclimation is also considered, that is the ability of microalgae to adjust their pigment content and composition under varying light and nutrient conditions. Photoacclimation was first represented as a steady state process, then the model was later extended to represent photoacclimation dynamics [9]. Likewise, Nikolaou et al. [10] have recently proposed a model building on both Han and Droop models to describe the dynamics of photoproduction, photoinhibition and photoacclimation. Multiphysics models that integrate these biological state models within computational fluid dynamics simulation in order to represent the effect of light attenuation and mass-transfer limitation, have also started to emerge [11, 12]. Dynamic models for long-term algae bioprocess simulation in both laboratory scale system [13] and industrial-scale system [14], which incorporate more operating factors (light attenuation, nutrient concentration, temperature, etc.), can also be found in the literature.

Recently, Nikolaou et al. [15] have proposed a semi-mechanistic model capable of quantitative prediction of the flux of fluorescence by taking into account three distinct processes acting on different time scales: photoproduction encompassing all the processes that are responsible for the capture and utilization of photons; photoinhibition, the observed loss of photosynthetic productivity as a result of excess or prolonged exposure to light that is associated with the damage of functional components of the photosynthetic apparatus; and photoregulation, also referred to as non-photochemical-quenching (NPQ), encompassing all the regulatory mechanisms that dissipate excess excitation energy as heat. It was established through this work that the dynamics of chlorophyll-a fluorescence can help provide good predictions of the photosynthetic response under variable light conditions, thus allowing for the mathematical modelling of key photosynthetic mechanisms. However, it has recently been found that this fluorescence model may not describe photoregulation adequately over long-term experiments, thus advocating for a more detailed and biologically consistent representation.

The development of reliable mathematical models presents many challenges, particularly when multiple phenomena spanning several orders of magnitude both in time and space scales are to be accounted for, as is the case in [15]. This results in a number of identifiability (and therefore experimental) issues [16]. For the purpose of design, control and optimisation nonetheless, model identification is key to attributing a biological meaning to the parameters as well as to enabling accurate predictions under conditions that differ from those used during the identification experiments (extrapolation) [17].

Experiments are often designed and planned based on the available biological insight about the system, but as the complexity of phenomena and their correlation (as typically occurs in multiscale systems) increases, an experimenter’s expertise alone may no longer be sufficient to identify a model in a statistically meaningful way. Accordingly, more informative experimental protocols are required, as determined for instance with an optimal design of experiments approach [18]. In particular, model-based design of experiments (MBDoE) provides a systematic approach to determining an experimental protocol that would maximise the level of information available for identifying that model [18]. This approach has been successful in numerous applications, including biological systems [19, 20].

The core of this contribution is to demonstrate how our understanding and prediction capability of photosynthetic process may benefit from a synergy between modelling and experiments using MBDoE. As a case study we consider the model by [15] as a starting point. Experimental evidence [21] suggests that photoregulation mechanisms involve at least two interdependent processes, which calls for a more complex NPQ representation in this model. At the same time, a more complex NPQ representation augments the number of unknown parameters, a condition that we will show can be compensated by a wider array of experimental measurements combined with MBDoE. Eventually, we will be able to propose a biologically consistent description of the NPQ mechanism leading to a significantly improved prediction capability of photosynthetic production under a great variety of light conditions. Finally, we will be able to verify the link between fluorescence and photosynthesis-irradiance (PI) curves experimentally, a link that was theoricized in [22].

The remainder of the article is organised as follows. Section 2 presents the experimental methods and protocols used to develop and identify the proposed model. The dynamic model of chlorophyll fluorescence is presented and discussed in Section 3.1, with a particular focus on the description of NPQ activation/relaxation. This new model calls for a tailored calibration procedure in order to distinguish between the dynamics of NPQ and photoinhibition, whose time-scales are overlapping. The model’s identifiability is analysed in Section 3.2. Section 4 applies this methodology to accurately calibrate and validate the model for the microalga N. gaditana. A first calibration of the model parameters is conducted in Section 4.1, followed by a sensitivity analysis and model-based design of experiment in Section 4.2 and 4.3, respectively. The calibration results are discussed in Section 4.4, and the model is validated against multiple additional experiments in Section 4.5. Finally, Section 5 concludes the paper. Symbols and acronyms used in this work are listed and defined in Table 1.

Table 1

Symbols and acronyms.

meaning.

Symbol	Meaning and units
αF	activity level of fast energy dependent quenching [-]
αS	activity level of slow energy dependent quenching [-]
αSS	reference function for energy dependent quenching activity [-]
A	fraction of photosynthetic units in open state
B	fraction of photosynthetic units in closed state
C	fraction of photosynthetic units in inhibited state
ηD	rate of basal termal decay relative to the rate of fluorescence [-]
ηI	rate of inhibition related quenching relative to the rate of fluorescence [-]
ηP	rate of photoproduction relative to the rate of fluorescence [-]
ηqE	rate of energy dependent quenching relative to the rate of fluorescence [-]
η¯qEC	maximum rate of interaction energy dependent quenching relative to the rate of fluorescence [-]
η¯qEF	maximum rate of fast energy dependent quenching relative to the rate of fluorescence [-]
η¯qES	maximum rate of slow energy dependent quenching relative to the rate of fluorescence [-]
F′	maximum rate of slow energy dependent quenching relative to the rate of fluorescence [-]
F0	dark-adapted minimal fluorescence flux
F0′	light-adapted minimum fluorescence flux
Fm	dark-adapted maximal fluorescence flux
Fm′	light-adapted maximum fluorescence flux
ΦfA	fluorescence quantum yield of a reaction centre in state A [-]
ΦPA	quantum yield of photosynthesis of an open reaction centre of the photosystem II [μmol_e^-μE^-1]
ΦfB	fluorescence quantum yield of a reaction centre in state B [-]
ΦfC	fluorescence quantum yield of a reaction centre in state C [-]
ΦfC	fluorescence quantum yield of a reaction centre in state C [-]
Φf	quantum yield of fluorescence [−]
ΦP	photosynthesis quantum yield [μmol_O_2μE^-1]
ΦPS2	realised quantum yield of photosynthesis [μmol_e^-μE^-1]
IqE	irradiance level at which half of the maximal qE activity is realised [μEm−2s−1]
kd	damage rate constant [-]
kr	repair rate constant [s−1]
N	chlorophyll specific number of photosynthetic units [μmol_O_2gchl-1]
n	Hill parameter related to the shape of sigmoid function describing NPQ activity [-]
ν	number of electrons per molecule of dissociated H2O
P	photosynthesis rate [μmol_O_2gchl-1s^-1]
qNPQ	NPQ index [-]
SF	scaling factor for fluorescence model (proportional to the chlorophyll content)
σ	total cross section [m^2gchl-1]
σPS2	effective cross section of photosystem II [m2 μE−1]
τ	turn over rate [s]
ξF	time constant of the fast NPQ activation/relaxation mechanism [s−1]
ξS	time constant of the slow NPQ activation/relaxation mechanism [s−1]
Acronym	Meaning
ASII	functional antenna size
DCMU	3-(3,4-dichlorophenyl)-1,1-dimethylurea
MBDoE	model-based design of experiments
NPQ	non-photochemical-quenching
PAM	pulse amplitude modulation
PI	photosynthesis-irradiance
PSII	photosystems II
PSU	photosynthetic unit
RCII	reaction centre of photosystem II

2 Material and Methods

2.1 Strain Cultivation

Microalgae Nannochloropsis gaditana (CCAP, strain 849/5) were grown in a sterile, filtered F/2 medium, using sea salts (32 gL−1) from Sigma, 40 mMTris HCl, pH 8 and Sigma Guillard’s (F/2) marine water enrichment solution. Growth experiments were performed in the multi-cultivator MC 1000-OD system (Photon Systems Instruments, Czech Republic) at a temperature of 21°C and a light intensity of 100 μEm−2s−1 provided continuously by an array of white LEDs. The suspension culture was constantly mixed and aerated by bubbling air. Pre-cultures were grown at 100 μEm−2s−1 in glass bottles of 0.25 L under a continuous airflow, enriched with 5% CO2. At the exponential phase, the pre-culture was centrifuged and re-suspended in fresh medium to reach a final concentration of 9 · 106cellmL−1, before introduction in the multi-cultivator.

Three types of measurements are used for model validation and calibration purposes later on, which are described in the following subsections.

2.2 PAM fluorometry

When a sample containing microalgae is exposed to light, a fraction of the incoming photons is absorbed by pigment molecules, another fraction is scattered out and the rest passes through. The absorbed photons have three possible fates: they can either be captured by the reaction centre of photosystem II (RCII) to drive photosynthesis (photoproduction), or dissipated as heat (photoregulation), or re-emitted as fluorescence [23]. Thus, much information about the photosynthetic processes can be inferred by measuring the fluorescence fluxes under specific lighting protocols that preferentially activate or inactivate the photoproduction and photoregulation mechanisms.

Today’s state-of-the-art equipment, such as pulse amplitude modulation (PAM) fluorometers, can implement complex light protocols and record fluorescence fluxes both fast and accurately. PAM fluorometry excites microalgae photosynthetic apparatus by using three distinct light sources: a weak modulated light to measure the fluorescence flux; an actinic light to drive photosynthesis by exciting the photosynthetic apparatus and activating photoregulation; and a light pulse of high intensity to saturate the photosystems II (PSII).

Depending on the light conditions, the following fluorescence fluxes can be recorded by a PAM fluorometer [24]:

Dark-adapted, minimum fluorescence flux, F0: The measuring light is applied to a sample that has been kept in the dark a sufficient length of time to completely oxidize the RCII. In the dark adapted sample the active RCII are thus ready to accept incoming photons (open state) and the NPQ processes are inactive.

Dark-adapted, maximum fluorescence flux, Fm: The measuring light is applied, after a short and intense light pulse, to a dark adapted sample. The PSII are completely reduced after the intense light pulse, so the active RCII are all occupied and cannot accept further photons (closed state), while the NPQ processes remain inactive.

Light-adapted, minimum fluorescence flux, : The measuring light is applied to a sample that has been exposed to a constant actinic light, after switching-off of the actinic light. The PSII are completely oxidized and the NPQ processes are active due to the application of an actinic light.

Light-adapted, maximum fluorescence flux, : The measuring light is applied, after a short and intense light pulse, to a sample that has been exposed to a constant actinic light. The PSII are completely reduced and the NPQ processes are active.

Light-adapted, steady-state fluorescence flux, F′: The measuring light is applied to a sample that has been exposed to a constant actinic light. The PSII are only partially reduced, meaning that the active RCII are in a mix of open and closed states and the NPQ processes are active.

All of the above-mentioned fluorescence fluxes also depend on the presence of damaged (inactive) RCII that result from the evolution of photoinhibition when the incoming photons flux exceeds the capacity of the PSII.

In this work, all the fluorescence measurements are performed using a Dual PAM (Walz, Germany), after a dark adaptation period of 20 min. For model calibration purposes three fluorescence experiments are performed. The first experiment (Exp1) considers a light protocol with increasing actinic light intensities applied in steps of 60 s. Before switching-on of the actinic light and during the final 2 s of each step, a saturating light pulse of 6000 μEm−2s−1 is applied during 0.6 s, followed by a dark period (actinic light off) of 1.4 s; measurements are recorded before and after the saturating pulses and after the dark periods, which correspond to F′, and respectively (see section 4.1). The second experiment (Exp2) considers an actinic light profile at 2000 μEm−2s−1 during 10 minutes followed by a dark period of 60 minutes. Saturating pulses are applied at regular intervals during this experiment in order to measure F′, and . The third experiment (Exp3) is optimally designed in order to improve upon the parameter estimation as described later in Section 4.4.

Our previous work [15] assumed a constant 1% variance of the fluorescence flux measurements. Here, we assume a variance model whereby a constant term (equal to 0.0015) is added to the 1% of the fluorescence fluxes, accounting for the sensitivity of the photomultiplier measuring the fluorescence signal. The utilization of a more realistic variance model for the measuring device may cause further practical identifiability problems.

2.3 Photosynthesis Rate

Out of the possible ways of measuring the photosynthesis rate, we consider measurements of the maximal rate of photosynthetic oxygen evolution at a specific actinic light using a Clark electrode (Hansatech, UK). In order to use PI measurements for parameter estimation, knowledge of the exact experimental light protocol is paramount, as established in Bernardi et al. [22]. Here, we consider PI measurements obtained from three separate experiments. Each experiment, involves exposing a dark adapted sample to two different light intensities over variable time periods (the exact protocols are reported in Table 2). The photosynthesis rate is measured at the end of each time period, thus providing a total of six experimental PI points. Sample 1 is to be included in the calibration data set, whereas the other two experiments are used for model validation only.

Table 2

PI experimental protocol.

Light irradiance and corresponding time duration for the PI measurements in all three experiments.

	Sample 1		Sample 2		Sample 3
irradiance [μEm−2s−1]	400	1500	100	750	250	3600
stage duration [s]	130	150	230	200	150	130
photosynthesis rate [g_o_2gchl-1h^-1]	4.10	5.58	1.22	5.18	2.58	7.22

2.4 Antenna Size Measurements

Antenna size measurements are used to study the saturation dynamics of the PSII. A LED pump and probe JTS10 fluorometer in fluorescence mode is used for the measurement. Fluorescence inductions are measured in the infrared region of the spectrum upon excitation with blue light at 450 nm. 3-(3,4-dichlorophenyl)-1,1-dimethylurea (DCMU) is added at a concentration of 80 μM to prevent oxidation of the primary quinone acceptor QA [25, 26]. In the presence of this inhibitor, the half-saturation time constant of the fluorescence rise response is known to be inversely proportional to the so called functional antenna size (ASII) [27].

We consider antenna size experiments for five different actinic light intensities (45, 80, 150, 320, 940 and 2050 μEm−2s−1) at 630 nm and four replicates are measured at each light intensity. Fig 1 shows a set of such experimental saturating curves. For each saturating curve the value of ASII is derived by fitting the data with a first order model, as further described in Section 3.2.

Fig 1

Antenna size experiment.

Evolution of the fluorescence flux of PSII from dark-adapted acclimated cells treated with DCMU. The time required for reaching half of the maximum response is inversely proportional to ASII. The fluorescence flux is reported for three different actinic lights and normalised by the maximum flux value (a.u.: arbitrary units).

2.5 Numerical Aspects

Simulations of the fluorescence model are conducted in the modelling environment gPROMS (Process System Enterprise, gPROMS v 4.1, www.psenterprise.com/gproms, 1997–2016), together with the formulation and solution of parameter estimation and MBDoE problems. Parameter estimation is performed using maximum likelihood estimation and statistical confidence analysis [28], in order for the model predictions to match the measured fluorescence fluxes (, and F′), and ASII and PI measurements. MBDoE optimises the actinic light profile and the measuring instants during a PAM experiment by minimizing the A-criterion [29], as other optimal criteria (namely D- and E-criteria) present numerical issues. The A-criterion aims to minimize the trace of the variance-covariance matrix of the parameter estimates to determine an information-rich experiment.

3 Model Development and Analysis

3.1 Modelling Chlorophyll Fluorescence

The model by Nikolaou et al. [15] predicts fluorescence fluxes by taking into account the state of the PSII and the activity of photoregulation. This model builds upon the well established state model by Han et al. [7] describing the state of the PSII in terms of open, closed and damaged RCII, and it introduces an activity equation describing the evolution of photoregulation in a semi-empirical manner (first-order process).

This latter modelling assumption introduces an important limitation since, from a biological point of view, NPQ in algae is related to at least two distinct processes. The first one is a fast process involving the LHCSR protein, with a time constant of seconds for both its activation and relaxation [30]. The second one is related to the xanthophyll cycle and acts on a time scale of minutes. More specifically, zeaxanthin can have a complex effect on the NPQ activity as it both enhances the quenching effect of LHCSR and acts as an additional quencher [21]. In the model by Nikolaou et al. [15] the two interdependent NPQ processes are lumped into a single first order process. This approximation of the actual biological mechanism may lead to poor prediction of the fluorescence fluxes, for instance in experiments where the actinic light is kept constant for several minutes. This is illustrated in Fig 2(a), where the model predictions obtained with the parameter values established in [15], significantly deviate from the data generated by an experiment with a constant actinic light applied for 600 s, followed by a recovery period of 3600 s. Moreover, Fig 2(b) suggests that the recovery phase involves a two time-scale NPQ relaxation mechanism rather than a simple first-order process.

Fig 2

Constant actinic light PAM experiment.

(a) Comparison between the predicted and measured fluorescence fluxes (triangles), (squares) and F′ (circles) in response to a constant light experiment. The grey-shaded area represents the light intensity. (b) Measured value of qNPQ, defined as , during the recovery phase of experiment Exp2 along with predicted values using different modelling assumptions. The dashed lines consider a first-order model to represent NPQ; the solid line considers the NPQ as the combined effect of two interdependent processes with different time scales.

Given that a first-order process is not sufficient to accurately describe the actual photoregulation mechanisms, we propose the following, more complex dynamics: The combined effect of LHCSR protein and zeaxanthin is apparent in Eq (4), where the parameters and represent the rates of the fast (LHCSR-related) and slow (zeaxanthin-related) NPQ processes, both relative to the rate of fluorescence. The parameter represents the enhancing effect of zeaxanthin in the quenching capability of the LHCSR protein. The activities of the fast and slow NPQ processes are described by the conceptual variables α and α, both modelled as first-order processes. Moreover, the reference α is the same in Eqs 1 and 2 since LHCSR- and zeaxanthin-related quenching are both triggered by low lumenal pH [21], yet with different time constant ξ and ξ. The reference α itself is modelled as a sigmoid (Hill) function, in agreement with experimental measurements of the NPQ index as function of I in [31], with IqE and n representing the irradiance level at which half of the maximal NPQ activity is triggered (α = 0.5) and the sharpness of the switch-like transition, respectively. Overall, the NPQ dynamic is parametrized by , , , ξ and ξ, instead of only two parameters in the original model [15]. The remaining model equations remain unchanged compared with [15], and are listed below for completeness:

Eq (5) defines the fluorescence flux, as a function of the total cross-section, ; the fluorescence quantum yield, Φ [−]; and a parameter depending on the characteristics of the PAM fluorometer and the chlorophyll content of the sample, S [Vgchlm−2]. The fluorescence quantum yield, given in Eq (6), is the harmonic mean of the fluorescence quantum yields of open (A), closed (B) and damaged RCII (C). The dynamics of A, B and C are described according to the Han model [7] in Eqs (7)–(9), with the following parameters: the effective cross-section of the PSII, σPS2 [m2 μE−1]; the turnover rate, τ [s−1]; the damage rate constant, k[−]; and the repair rate constant k[s−1]. Eq (10) relates the effective cross-section, σPS2 of the Han model with σ, effectively linking the state of the RCII with the fluorescence flux. It is important to note that while σPS2 is a constant in the Han model, it becomes a function of the photoregulation activity (Eq 4), via the quantum yield of photosynthesis of an open RCII (Eqs 11 and 12). The parameter N represent the chlorophyll specific number of photosynthetic units , and ν is a stoichiometric factor reflecting the minimum theoretical value of 4 electrons for each dissociated water molecule per the water dissociation reaction (2H2O + 4e− → O2 + 4H+). Eq (11) defines the quantum yield of photosynthesis of an open RCII. Finally, Eqs (12)–(14), define the quantum yield of fluorescence of the open, closed and inhibited PSUs, with the parameters η, η and η representing, respectively, the rates of photoproduction, basal thermal decay in dark-adapted state and qI-quenching respectively. Overall, the maximum and minimum fluorescence fluxes and can be calculated from Eqs (5) and (6) by varying A and B. Specifically is obtained by imposing A = 0 and B = 1 − C, whereas by imposing A = 1 − C and B = 0. Moreover, the distinction between dark and light-adapted fluxes depends on the value of the variables α and α, with α = α = 0 for the dark-adapted state.

The improved representation of the NPQ mechanisms in Eqs (1)–(4) increases the complexity of the model in terms of the number of parameters from 13 to 16. The augmented parameter vector is as follows: The parameter estimation problem is further complicated by the fact that the three recovery processes are now overlapping, while NPQ relaxation and damaged PSUs recovery were acting on distinct time scales in the original model [15].

3.2 Coupling Fluorescence with PI and ASII Measurements

The models (1)–(14) presents some identifiability issues [32] in the sense that not all of its parameters can be estimated uniquely when only fluorescence measurements are available. We note that, even in the original model [15], the parameters η, N, k and τ could not be uniquely estimated.

Based on the work by Falkowski and Raven [33], Bernardi et al. [22] derived an expression of the photosynthesis rate as: with Φ the photosynthesis quantum yield []. The units for p are [], which by definition is the chlorophyll-specific photosynthesis rate expressed in terms of O2 production. The fluorescence models (1)–(14), predicts the value of the realized quantum yield of photosynthesis, in units of [], from which we then derive [34]. This way, the photosynthesis rate expression (15) does not introduce any additional parameters and conveniently couples PI measurements with fluorescence measurements in order to conduct parameter estimation. Although parameters σ, S and N are structurally unidentifiable when using fluorescence measurements only, they become identifiable in adding PI measurements to the calibration set.

Note that fluorescence or PI measurements alone may only give information about the steady-state values of the fraction of open/closed RCII as a function of the light intensity, since the measurement frequency is not high enough to probe their dynamics. Based on such steady-state information it is only possible to estimate the ratio between σPS2 and τ. In contrast, ASII measurements can provide information about this fast time scale, thereby providing a means of estimating the parameters σPS2 and τ independently.

The addition of DCMU to a dark-adapted sample inhibits the transition B → A, by preventing oxidation of the primary quinone acceptor Q. This way, the fluorescence induction curves shown in Fig 1 contain information about the dynamics A → B only. The time constant of this activation process can be obtained by fitting a first-order response model [6-8] for each light intensity.

The values of ASII for the different light intensities are reported in Fig 3 along with a regression fit. The standard deviations remain small with a maximum value of about 2.5% of the ASII value at 2050 μEm−2s−1. The slope of the regression line is 0.25 ± 0.009[μE−1m2] with an R2 value of 0.993.

Fig 3

Antenna size experiment.

Linear regression of ASII from the antenna size measurements at five different light intensities. The shaded area represents the confidence region of the linear regression.

Fig 3 validates the Han hypothesis of a linear correlation between ASII and the light intensity, and provides useful information to enhance the parameter estimation of the fluorescence model. Specifically, in our models (1)–(15) the antenna size corresponds to the product between the effective cross-section σPS2 and the light intensity I, ASII = σPS2 I. Therefore, the slope of the linear regression in Fig 3 represents the parameter of the Han model σPS2 for a dark adapted sample (α = α = 0). Recall also that in the models (1)–(15), σPS2 is related to the parameter σ via Eq (10) and is a function of the photoregulation activity. In particular, the time horizon of 0.35 seconds for the antenna size experiment is sufficiently short to prevent significant activation of the photoregulation mechanisms. One way to exploit the information obtained through the antenna size experiments consist of adding a virtual experiment to the calibration set, whereby the value of σPS2 is a measured variable with a standard deviation of 0.009 μE−1 m2.

4 Model Validation Methodology

4.1 Step 1: Preliminary Identifiability Analysis

The model parameters are estimated using the two PAM experiments Exp1 and Exp2 described in Section 2.2, the PI measurements from Sample 1 described in Section 2.3 and the ASII measurements in Section 2.4. It can be established, e.g. by using a differential algebra approach as available through the software package DAISY [35], that all the parameters but η are structurally identifiable with such a calibration set. In principle, one could estimate the parameter η, representing the rate of basal thermal decay relative to the rate of fluorescence, based on the probability of thermal dissipation and the probability of fluorescence for a photon absorbed by a dark adapted RCII. Because this kind of measurements was not available to us, we set this parameter to η = 5 subsequently, which is the mean of those η values for which the computed fluorescence quantum yields are consistent with the data from Huot and Babin [36].

The results reported in Table 3 show that part of the estimated parameter values are not statistically meaningful. In particular:

Table 3

Preliminary model calibration.

Parameter estimates along with their 95% confidence interval and t-values. The reference t-value is 1.65. The calibration set is comprised of Exp1, Exp2, ASII measurement and Sample 1 of PI measurements.

Parameter	Estimated value	95% conf. int.	t-value 95%	Units
ξF	1.83 × 10−1	1.91 × 101	0.0096 *	s−1
ξS	9.68 × 10−4	6.57 × 10−5	14.74	s−1
IqE	5.96 × 102	4.18 × 101	14.26	μEm−2s−1
kd	2.04 × 10−6	1.32 × 10−6	1.55 *	−
kr	2.78 × 10−5	5.98 × 10−5	0.46 *	s−1
N	5.31 × 10−1	8.61 × 10−2	6.16	μmol_O_2gchl-1
n	2.18 × 100	1.74 × 10−1	12.54	−
ηI	2.77 × 100	1.45 × 100	1.92	−
η¯qEF	8.17 × 100	9.18 × 10−1	8.90	−
η¯qES	1.92 × 101	1.20 × 100	16.02	−
η¯qEC	2.44 × 101	3.32 × 100	7.35	−
ηP	1.14 × 101	3.07 × 10−1	36.97	−
σ	7.79 × 10−1	1.19 × 10−1	6.55	m^2gchl-1
τ	8.45 × 10−3	1.07 × 10−3	7.87	s
SF1 a	1.77 × 100	2.72 × 10−1	6.53	Vgchlm−2
SF2 a	1.97 × 100	3.04 × 10−1	6.47	Vgchlm−2

a S refers to Exp1, S refers to Exp2 experiment. The different values are due to different cell concentrations in the respective samples.

* an individual 95% t-value smaller than the reference t-value indicates that the available data may not be sufficient to estimate the parameter precisely

the parameters k and k may not be confidently estimated due to their large correlation with η;

the parameter has a very low sensitivity caused by the fact that the dynamics of LHCSR protein activation are much faster than the resolution of the PAM fluorometer.

It is well known that even when a model is structurally identifiable, can measurements noise and other sources of uncertainty still hinder its practical identifiability [19, 37]. Since the parameters k, k, and η are highly correlated, it is unlikely that generic experiments will provide the level of information that is needed to achieve statistically significant estimates for these parameters.

4.2 Step 2: Sensitivity Analysis

The sensitivity of a parameter θ with respect to a measured variable y can be computed as (y′ − y)/Δθ, where y and y′ are the measured variables calculated with the nominal and perturbed value of the parameter, respectively, and Δθ is the parameter perturbation. By calculating the sensitivity of each parameter one can derive useful information about a model’s practical identifiability.

We illustrate this technique by analysing the sensitivity of the parameter ξ with respect to the measurement times of the fluorescence fluxes, considering a 5% increase in the value of ξ. We simulate an experiment in which the actinic light is switched on for 60 s at 2000 μEm−2s−1 and then switched off for another 60 s. The time-varying sensitivity profiles of , and F′ with respect to the parameter ξ are plotted in Fig 4. The sensitivities of the three fluxes pass by a maximum within a few seconds after the light switch occurs and then decrease exponentially to zero, confirming that ξ may only be estimated if a measurement is taken shortly after the light has been switched on or off.

Fig 4

Sensitivity analysis.

Sensitivity trajectories, s, of the fluxes , and F′ with respect to the parameter ξ. The red continuous line represents the sensitivity of , the green dotted line represents the sensitivity of F′, and the blue dashed line represents the sensitivity of . The light protocol used to obtain these curves is 60 seconds of actinic light at 2000 μEm−2s−1 followed by 60 seconds of dark. The nominal value of ξ is 0.18 s−1.

In a standard PAM experiment, like the two experiments used for calibration in the previous section, measurements are usually taken within 40 to 60 s after a variation in actinic light. It is therefore clear that the sensitivity of the measurements will usually be very low with respect to the parameter ξ and the resulting estimate may not be confident. Instead, an informative (designed) experiment should ensure that some of the measurements are sampled when sensitivities of the measured variables are indeed maximal.

4.3 Step 3: Model-Based Design of Experiments

MBDoE is an effective methodology for designing experiments that contain maximal information based on a mathematical model, whose parameters are to be estimated. The experiment design is formulated as an optimisation problem wherein the decision variables correspond to time-invariant or time-varying inputs and the measurement times for the variables [18].

In practice, experimental designs often fail to take into account the relationship between the measured variables and the model parameters in order to effectively excite the system and to collect the measurements where they are the most informative. In other words, experimental protocols for model identification are not always designed in order to guarantee that the vector y(t) of measured output variables is maximally sensitive to the model parameters, as could be observed in the previous section for the parameter ξ.

Here, the measured outputs are the three fluorescence fluxes , and F′ that can be optimally sampled at 57 points during the experiment, with two consecutive measurements being separated by at least 40 s, a constraint posed by biological considerations. The experimental design variables consist of the corresponding measurement times, along with the actinic light profile as discretized on 20 constant light stages and whose durations are also decision variables. Note that the PAM fluorometer used to conduct the experiments can only apply discrete actinic light intensities of 0, 6, 13, 22, 37, 53, 70, 95, 126, 166, 216, 273, 339, 430, 531, 660, 825, 1028, 1287, 1594 and 1952 μEm−2s−1. The total experimental horizon has been set equal to 2400 s and at the end of the experiment a dark period of 3 hours is added to better probe the slow dynamics of damaged RCII recovery. To summarise, the design vector includes: the measuring instants of the three fluorescence fluxes, the durations of actinic light stages, and the level of actinic light during each stage.

The new fluorescence experiment designed with MBDoE, based on the preliminary parameters values obtained in Section 4.1 is shown in Fig 5(c), and referred to as Exp3 subsequently. Notice how MBDoE leads to a complex dynamic light protocol so as to excite the system in order to provide a highly informative input-output data set.

Fig 5

Final calibration results.

PAM fluorescence profiles along with the model predictions corresponding to Exp1, Exp2 and Exp3 for the model parameters in Table 4.

Table 4

Final parameter estimation.

Parameter estimates along with their 95% confidence interval and t-values. The reference t-value is 1.65. The calibration set is comprised of Exp1, Exp2, Exp3, Sample 1 of PI measurements and ASII measurements.

Parameter	Estimated value	95% conf. int.	t-value 95%	Units
ξF	2.68 × 10−1	3.50 × 10−2	7.67	s−1
ξS	1.32 × 10−3	6.97 × 10−5	18.88	s−1
IqE	5.95 × 10 2	2.07×101	28.76	μEm−2s−1
kd	9.95 × 10−7	2.67 × 10−7	3.73	−
kr	5.10 × 10−5	2.67 × 10−5	1.78	s−1
N	4.83 × 10−1	7.52 × 10−2	6.43	μmol_O_2gchl-1
n	2.40 × 100	1.27 × 10−1	18.87	−
ηI	1.41 × 101	3.98 × 100	3.54	−
η¯qEF	5.96 × 100	4.98 × 10−1	11.95	−
η¯qES	1.23 × 101	5.75 × 10−1	21.35	−
η¯qEC	2.47 × 101	1.69 × 100	14.58	−
ηP	1.04 × 101	2.33 × 10−1	44.54	−
σ	7.33 × 10−1	7.50 × 10−2	6.84	m^2gchl-1
τ	6.95 × 10−3	7.50 × 10−4	9.26	s
SF1 a	1.81 × 100	3.01 × 10−1	6.82	Vgchlm−2
SF2 a	2.06 × 100	3.01 × 10−1	6.81	Vgchlm−2
SF3 a	1.30 × 100	1.90 × 10−1	6.82	Vgchlm−2

a S refers to Exp1; S refers to Exp2; S refers to Exp3 experiment. The different values are due to different cell concentrations in the respective samples.

4.4 Step 4: Final Parameter Estimation

The designed experiment (Exp3) is included along with Exp1 and Exp2 in the parameter estimation problem. The updated parameter estimates along with their 95% confidence intervals and t-values are listed in Table 4, and the fitting results of the model are presented in Figs 5 and 6.

Fig 6

Experimental PI measurements along with model prediction.

The solid line represents the model predicted PI curve; the black squares are the experimental data used for model calibration; and the blue stars represent the experimental data used for model validation.

The addition of experiment Exp3 now allows for the accurate estimation of both ξ and k, and besides, k is also estimated in a statistically meaningful way. The rest of the parameters are also obtained with better confidence, thus demonstrating the effectiveness of MBDoE techniques. The estimated value of η is significantly different from the one obtained in Table 3, which is not surprising considering its high correlation with the parameters k and k, and the fact that such parameters could not be properly estimated in the preliminary calibration in Section 4.1.

As shown in Fig 5, the model capability to represent the experimental data is generally very good, although some small mismatch can still be noticed. The first mismatch is related to the recovery phase of the three experiments, as the value of is underestimated by the model in Fig 5(a) and 5(c), whereas it is overestimated in Fig 5(b). Here, we have to consider that the three experiments were conducted using samples drawn from different cultures. Although the growth conditions were identical for all three cultures, the intrinsic biological variability of the system is likely to be responsible for the different responses in the experimental data.

The second mismatch is evident in Fig 5(c) between 1000 and 2000 s, where the predicted flux is underestimated during the low light periods. Although it is systematic this underestimation remains quite small. Future experiments will help understand whether this mismatch is related to structural inconsistency of the model, or if it can be attributed to biological variability or experimental issues.

The PI measurements simulated by the model along with the experimental values are reported in Fig 6. The black squares represent the data used for model calibration, whereas the blue stars are considered for model validation later in Section 4.5. The model shows a good agreement with the measured PI values used for model calibration (Sample 1). An interesting finding is that only two PI measurements are sufficient here to identify the parameter N, thus removing the need for time-demanding PI measurements.

Although some small discrepancies are observed, which may relate to the intrinsic variability of the microalgae sample, the model is able to reproduce experiments that are very different from each other, and statistically meaningful parameter estimates can be determined for all but one (η) parameters.

4.5 Step 5: Model Validation

Three additional PAM experiments are considered in order to validate the model. The validation set includes one constant light experiment (Val1) and two different variable light experiments (Val2 and Val3) where the sample is subject to complex actinic light profiles. The predicted fluorescence fluxes with the parameter estimates given in Table 4 along with the experimental data are shown in Fig 7. Moreover, two additional PI experiments are used for validation purposes, each consisting of two PI measurements, as shown in Fig 6.

Fig 7

Model validation results.

PAM fluorescence profiles along with the model predictions corresponding to Val1, Val2 and Val3, for the model parameters in Table 4.

The model predicts all three validation experiments in a very satisfactory way and is capable of capturing the variations in photosynthetic response triggered by different light dynamics. Moreover, the model gives accurate prediction of photosynthesis rate measurements. In particular, it is important to note that the model predictions remain accurate even for PAM protocols very different from the ones used in the calibration set, that is when the model is used for extrapolation. This confirms that the model provides an effective tool to predict PI curves without the use of classical measuring techniques that take time and resources to implement. The main discrepancy between the experiments and the model predictions can be observed in Fig 7(c) during the first 300 s of the experiment, which is similar in nature to the mismatch previously observed in Fig 5(c). However, the fits are generally good, thus confirming that the model is capable of accurate quantitative prediction for a wide range of PAM experiments. To our knowledge, no other mathematical model to date has been tested through such challenging experimental protocols or has demonstrated such a consistent prediction capability.

5 Conclusions

A mathematical model incorporating a detailed and biologically consistent representation of the NPQ mechanism has been proposed and validated using experimental data. The challenge of model identification has been tackled by introducing specific measurements, namely PI and antenna size measurements and by designing a tailored PAM protocol using MBDoE techniques. The results show that the model is able to predict chlorophyll fluorescence and photosynthesis rate with a good accuracy under a large variety of light conditions for N. Gaditana, thus paving the way towards a more reliable and realistic description of the effects of light dynamics on microalgae growth. Furthermore, the connection between PI curves and fluorescence has been verified practically, which could prove beneficial for reducing the experimental effort relative to PI curves experiments.

Future work will aim at moving up in the multiscale description of photosynthetic mechanisms and microalgae growth rate. This will involve incorporating a suitable description of photoacclimation phenomena into the model, following the work by Nikolaou et al. [10]. The focus will also be on combining the model with physics-based hydrodynamic and light attenuation models for full-scale photobioreactor simulation and optimisation [2].

Experimental data.

Excel file containing all the experimental data used in the paper.

(XLSX)

Click here for additional data file.