Enhancing control systems of higher plant culture chambers via multilevel structural mechanistic modelling

Modelling higher plant growth is of strategic interest for modern agriculture as well as for the development of bioregenerative life support systems for space applications, where crop growth is expected to play an essential role. The capability of constraint-based metabolic models to cope the diel dynamics of plants growth is integrated into a multilevel modelling approach including mass and energy transfer and enzyme kinetics. Lactuca sativa is used as an exemplary crop to validate, with experimental data, the approach presented as well as to design a novel model-based predictive control strategy embedding metabolic information. The proposed modelling strategy predicts with high accuracy the dynamics of gas exchange and the distribution of fluxes in the metabolic network whereas the control architecture presented can be useful to manage higher plants chambers and open new ways of merging metabolome and control algorithms.

1 Introduction

Modelling crop growth has been a topic of research since the mid-twentieth century given the relevance that agronomic-related activities have in the global economy, but the focus on plant modelling research has evolved in the last decades moving towards new applications (Louarn and Song, 2020). A lot of attention has been placed on developing full-canopy models to assess climate change from different perspectives, such as its effect on crop physiology, the higher plant adaptive strategies, or the contribution of forestry and agricultural systems on carbon dioxide capture, to mention a few (Soussana et al., 2010; Peng et al., 2020). Besides responding to the current global climate and demographic challenges, the need for more efficient forms of horticulture to increase productivities, improve yield, and optimize crop growth has also contributed to the generation of mathematical models to support agronomic activities. Still, much progress is required in the field of biological systems modelling and this is especially relevant in the case of higher plants due to the complexity of their underlying growth mechanisms. Modelling complex systems like higher plants may be an objective by itself as before-mentioned, but they may also lead to other interesting applications like the development of model-based control methodologies. Particularly, space research and the development of bio-regenerative life support systems (BLSSs), which are the set of technologies designed to guarantee life in long-term crewed missions (Eckart, 1995), have exploited the use of controllers based on first principles in opposition to surrogated and reduced order models, to improve both the management of missions under operation as well as to improve the design of future missions (Fulget et al., 1999). One of the main actors on BLSS research is the Micro-Ecological Life Support System Alternative (MELiSSA), a European Space agency (ESA) BLSS program, which is devoted to developing technologies accompanied by a modelling framework to support their research and development activities. One remaining task, which is in turn one of the most relevant consensuses of the development of BLSSs shared by the major space agencies, is the importance of developing technologies and mathematical models to grow plants on space, which are expected to be the major source of edible biomass in BLSSs (Gitelson and Lisovsky, 2002; Poughon et al., 2009; Dong et al., 2017).

Even though current models can cope with the evolution of biomass, the compounds involved in photosynthesis and respiration (O₂ and CO₂), and the nutrient uptake by the roots, most of the phenomena involved in plant growth are not addressed given their complexity and lack of knowledge. This complexity is associated with several factors; higher plants are multicellular, compartmentalized organisms undergoing strong metabolic changes associated with the cyclic switch between light and dark phases of the day. They are also characterised by having a complex substrate partitioning strategy with different parts being coordinated to uptake and distribute specific compounds. This has contributed to the preference of empirical at the expense of mechanistic models due to the usually satisfactory information provided by the former, especially in nominal conditions (Boscheri et al., 2012; Amitrano et al., 2020). Empirical models cover a limited range of operating conditions though; thus, the scope of their use is narrow and cannot contemplate all scenarios that plant culture may undergo. As an alternative, mechanistic approaches have also been deployed to understand the first principles behind higher plant growth, surpassing the Farquhar model (Farquhar et al., 1980) to calculate photosynthesis rates based on the enzyme kinetics of principal metabolic pathways. However, out of the photosynthesis process, many key mechanisms like respiration, substrate accumulation and management, tissue morphology, or multi-tissue interactions have not completely been mathematically characterised yet. To treat such a complex system, the use of metabolome information has attracted the attention of plant modellers as an alternative to gather multiple biological reactions formalized as a constraint-based metabolic model (Gomes et al., 2015). It should be highlighted that constraint-based metabolic models have been demonstrated to successfully address the plant diel cycle with a light phase with resource accumulation and a dark phase with resource depletion (Cheung et al., 2014), a critical phenomenon in higher plant metabolism very difficult to deal through first-principle approaches thus far. Several efforts have been recently placed on integrating available mechanistic information and omics data in a common multiscale modelling framework that could potentially be used by the plant computational biology research community to feed data in a single converging platform (Marshall-Colon et al., 2017; Xiao et al., 2017).

In this study, the modelling of higher plants is approached through the design of a multilevel organization of the mechanistic processes that take place during crop growth. To present the results, L. sativa has been used as an exemplary higher plant. The higher level in the hierarchy copes with the mechanistic phenomena corresponding to a higher characteristic length (i.e., crop chamber scale), whereas the lower the level the smaller the characteristic length of the modelled phenomena (i.e., enzyme rate). Information follows a top-bottom flow, and it is eventually used to calculate the metabolic flux distribution by applying a flux balance analysis (FBA). This multilevel modelling approach is firstly validated with experimental data and secondly integrated in a model-based predictive control, representing, to the best of the authors, the first attempt to incorporate cell metabolism in an advanced control strategy.

Overall, the modelling and control methodology presented in this study may pave the way for a more efficient and sustainable agriculture either for intensive cultivation systems or as a part of BLSSs in space exploration.

2 Model proposal

2.1 Multilevel mechanistic model

The model developed here is organized following a multilevel approach, considering the different levels of the plant, from canopy to metabolic level, and uses the output of the higher levels as the input to the lower levels. In this section, the models used in the different levels are explained. A graphical description of the model organization is presented in Figure 1 and detailed in the following paragraphs. For the sake of clarity, model parameters are described in Table 1 and constant parameters in Table 2.

FIGURE 1

Figure 1 Graphical scheme of the multilevel mechanistic model approach used for L. sativa prediction growth.

TABLE 1

Table 1 Model parameters.

TABLE 2

Table 2 Model constants.

2.1.1 Level 3: Modelling canopy growth

When dealing with whole-leaf or canopy modelling, it is necessary to consider the effect of shading among leaves, which is not accounted in single-leaf models (Poulet, 2018). The common practice is to use the leaf area index (LAI) as an indicative parameter of the leaf area density over ground surface. The photon flux density inside the canopy I declines along the canopy exponentially and is a function of the leaf area index:

$\begin{array}{ll}I=I_u{e}^{(-kLAI)}& (1)\end{array}$

Parameter I_u represents the photon flux density at the top of the canopy, and k represents the extinction coefficient. Extending (1), leaves can receive direct photon flux density (I_s) or diffuse photon flux density (I_d) as stated by Thornley (2002), which are expressed in terms of μmol m^-2 leaf s^-1 by using the extinction coefficient k:

$\begin{array}{ll}{I}_s=k·f_s·I_u+k·(1-f_s)·I_u·e^{(-k·LAI)}& (2)\end{array}$

$\begin{array}{ll}{I}_d=k·(1-f_s)·I_u{e}^{(-k·LAI)}& (3)\end{array}$

Notice in (2) that the parts of the canopy under direct irradiance also receive diffuse irradiance. The LAI term should be differentiated between the fractions exposed to direct and diffuse light sources (LAI_s and LAI_d respectively), as suggested by Thornley (2002):

$\begin{array}{ll}LA{I}_s=(1-e^{-kLAI})/k& (4)\end{array}$

$\begin{array}{ll}LA{I}_d=LAI-LA{I}_s& (5)\end{array}$

The derivative of LAI_s and LAI_d should be obtained and used to integrate I_s and I_d to calculate the overall irradiance received by the canopy:

$\begin{array}{ll}dLA{I}_s=e^{-kLAI} dLAI& (6)\end{array}$

$\begin{array}{ll}dLA{I}_d= ( 1-e^{-kLAI}) dLAI& (7)\end{array}$

$\begin{array}{ll}{I}_{s,l}=\underset{0}{\overset{LAI}{\int }}{I}_{s}dLA{I}_s& (8)\end{array}$

$\begin{array}{ll}{I}_{d,l}=\underset{0}{\overset{LAI}{\int }}{I}_{d}dLA{I}_d& (9)\end{array}$

The way light irradiates exposed and shadowed leaves strongly affects the energy balance in the leaf surface, with the shadowed leaves irradiated by diffused light mainly of a wavelength of 600 nm corresponding to the green-colour spectrum of transmitted light. Energy received by irradiance contains the direct (E_hs) and diffuse (E_hd) terms, which are expressed as follows:

$\begin{array}{ll}{E}_{h_s}=I_s·Na·c·h·LA{I}_s\sum _{i={\lambda }_{min}}^{i={\lambda }_{max}}\frac{{\gamma }_i}{{\lambda }_i}& (10)\end{array}$

$\begin{array}{ll}{E}_{h_d}=I_d·Na·c·h·LA{I}_d\frac{1}{{\lambda }_{600nm}}& (11)\end{array}$

Na represents the Avogadro number, c represents the velocity of light, and γ is the fraction of wavelength λ that compose the light directly irradiating the canopy. The total energy irradiated to the leaves is the summation of both Equations (10) and (11). The radiation energy emitted by the plants (E_r), the energy lost by convection (E_conv), and the energy lost by transpiration (E_tr) are determined by the following equations:

$\begin{array}{ll}{E}_r=ϵ\sigma \left(T_{leaf}^4-T_b^4\right)& (12)\end{array}$

$\begin{array}{ll}{E}_{conv}=C_p{k}_t\frac{P_b}{R{T}_b}(T_{leaf}-T_b)& (13)\end{array}$

$\begin{array}{ll}{E}_{tr}={\lambda }_{mol}E{x}_{H_{2}O}& (14)\end{array}$

In (12), ϵ and σ represent the leaf emissivity and the Stefan–Boltzmann constant, respectively. In (13), C_p and k_t represent the molar specific heat capacity at constant pressure and 298.15 K and the heat transfer coefficient, respectively, the latter being a function of the diffusion coefficient and the boundary layer thickness, as follows (see Supplementary S1 for details on D^t calculation and for an extended description of the boundary layer model):

$\begin{array}{ll}{k}_t=\frac{{\text{D}}^t}{\delta }& (15)\end{array}$

In (14), λ_mol is the water latent heat of vaporization and Ex_H2O the transpiration rate defined in the following section. Finally, T_leaf is the leaf surface temperature, T_b is the bulk temperature, and R is the ideal gas constant.

2.1.2 Level 2: Modelling gas exchange rates

This level is dedicated to calculating the uptake and release rates between the atmosphere and the leaves concerning the exchange gases (Ex_CO2 , Ex_O2 , and Ex_H2O ). The approach to modelling gas exchanges between the leaves and the atmosphere follows Fick’s law, with the concentration gradient being the driver of the molecular transport:

$\begin{array}{ll}E{x}_{C{O}_2}=\frac{{\text{G}}^{C{O}_2}}{P_b}(P_b^{C{O}_2}-P_l^{C{O}_2})& (16)\end{array}$

$\begin{array}{ll}E{x}_{O_2}=\frac{{\text{G}}^{O_2}}{P_b}(P_b^{O_2}-P_l^{O_2})& (17)\end{array}$

$\begin{array}{ll}E{x}_{H_{2}O}=\frac{{\text{G}}^{H_{2}O}}{P_b}(P_b^{H_{2}O}-P_l^{H_{2}O})& (18)\end{array}$

In (16)–(18), the atmospheric partial pressure is calculated assuming gases behave following the general gas equation, whereas the conductance for the different gases (G^CO₂ , G^O₂ , and G^H₂O ) and the internal (i.e., leaf) partial pressure ($P_l^{C{O}_2}$, $P_l^{O_2}$, $P_l^{H_{2}O}$) as well as the leaf area (LA) are calculated according to Poulet et al. (2018) and explained in Supplementary S1.

2.1.3 Level 1: Modelling enzyme kinetics

The gross photosynthesis rate (P_g) is calculated using the Farquhar model (Farquhar et al., 1980), which has been widely used to model photosynthesis phenomena (Harley et al., 1992; Morgan and Robles, 2002; Arnold and Nikoloski, 2011). In this model, P_g is determined through finding the limiting rate of the photosynthesis, which is caused by either the regeneration of ribulose 1,5-biphosphate (J), substrate of the ribulose-1,5-biphosphate carboxylase (RuBisCo), or the RuBisCO carboxylation rate itself (V_c).

On the one hand, the regeneration of ribulose 1,5-biphosphate depends on the potential rate of the light electron transport chain (F_LETC) and its capacity to generate reducing power. Thus, it is necessary to define an expression for F_LETC, which can be approximated by a quadratic equation (Nikolov et al., 1995):

$\begin{array}{ll}{F}_{LETC}=\frac{(F_{LETC}^{max}+\phi I_u)-\sqrt{{(F_{LETC}^{max}+\phi I_u)}^2-4\phi \theta I_u{F}_{LETC}^{max}}}{2\theta }& (19)\end{array}$

$\begin{array}{ll}\text{φ}=\frac{(1-\text{f})}{2}& (20)\end{array}$

In (19), the maximum F_LETC is represented by $F_{LETC}^{max}$ (see Supplementary S1 for calculation details). The photosynthetic photon flux density used is represented by I_u and θ is a convexity coefficient. The efficiency of energy conversion is represented by φ, which is a function of the fraction of absorbed photon flux unavailable for photosynthesis (f). As previously demonstrated, it is necessary to consider the direct and diffuse irradiances when considering the whole canopy. Therefore, combining with (4) and (5):

$\begin{array}{ll}{F}_{LETC(I_s,I_d)}=\underset{0}{\overset{LAI}{\int }}\frac{F_{LETC}(I_s) dLA{I}_s}{2}+\underset{0}{\overset{LAI}{\int }}\frac{F_{LETC}(I_d) dLA{I}_d}{2}& (21)\end{array}$

Considering that two electrons are necessary per molecule of NADPH generated, the light electron transport chain rate resulting from (21) is divided by 2. This is necessary to provide information to the FBA with consistent units considering the stoichiometry matrix used (see S2 for a full list of reactions and stoichiometry matrix).

As mentioned before, in Farquhar et al.’s (1980) approach, it is necessary to convert F_LETC defined in (21) to a flow of ribulose 1,5-biphosphate (Ru5P) regeneration (J), through the following expression derived from the stoichiometry of the light electron transport chain and the Calvin cycle (Farquhar et al., 1980):

$\begin{array}{ll}J=\frac{F_{LETC}}{2(2+2\varnothing )}& (22)\end{array}$

The RuBisCo carboxylation (V_c) kinetics is of Michaelis–Menten type and is a function of the leaf internal oxygen (O_i) and carbon dioxide (C_i) concentrations:

$\begin{array}{ll}{V}_c=V_{cmax}\frac{(C_i-\Gamma )}{C_i+K_c(1+\frac{O_i}{K_o})}& (23)\end{array}$

In (23), V_cmax represents the maximum carboxylation velocity of RuBisCo, K_c and K_o are the Michaelis–Menten half-saturation constants for the carboxylation and oxygenation activities of RuBisCo, respectively, and Г represents the carbon dioxide compensation point. Details about the calculation of the internal carbon and oxygen concentrations are found in Supplementary S1, based on the boundary layer approach defined by Poulet et al. (2018; 2020). V_c can be directly used to feed the FBA as an upper bound (see Figure 1).

Finally, the gross photosynthesis rate (P_g) is determined by finding the minimum between the Ru5P regeneration rate (J), the RuBisCo carboxylation rate (V_c), and the gas exchange rate (Ex_CO2 ) and the net photosynthesis rate (P_n) is determined by retrieving the RuBisCo oxygenation (V_o) to P_g as follows:

$\begin{array}{ll}{P}_g=min(V_c,E{x}_{C{O}_2},J)·LA& (24)\end{array}$

$\begin{array}{ll}{V}_o=V_c·\frac{V_{omax}}{V_{cmax}}·\frac{\frac{O_l}{K_o}}{\frac{C_l}{K_c}}& (25)\end{array}$

$\begin{array}{ll}{P}_n=P_g-V_o·LA& (26)\end{array}$

Notice that Equation (24) can be formalized due to the conversion applied in Equation (22) considering that J represents the regeneration rate of Ru5P and V_c represents its carboxylation rate. Parameter A indicates the surface of the crop growing area. In this study, the use of a metabolic matrix makes the conversion from light electron transport flux (F_LETC) to RuBP regeneration (J) unnecessary, because this information is already included in the stoichiometric matrix. Similarly, the discontinuity introduced by Farquhar et al. (1980) in (24) can be prevented with an FBA formulation as addressed in the following section.

2.1.4 Level 0: Stoichiometry matrix

In this level, a simplified network model of the photosynthetic leaves’ metabolism of L. sativa is described. The stoichiometry matrix is based on the work of Sasidharan (2012), which contains the distinguishing characteristics of L. sativa, such as the reduced starch content to store carbon and the definition of the elemental composition that makes up the macromolecules of the biomass. This model though has been extended to include relevant reactions like the pentose phosphate pathway or the photorespiration cycle, originally missing. In the model used, the cellular organelles are described as different compartments. The model also describes dark and light phases of the day by duplicating each one of the reactions. Hence, a diel model is achieved where both light and dark phases of the day account for separate pools of metabolites and organelle compartments. Here, only those metabolites that have been reported to be accumulated in one phase and consumed in the other are connected by exchange reactions that simulate the transference of nutrients between phases. For example, sugars that are synthesised in the light phase can be used in the dark phase due to the addition of exchange reactions among day phases. In Figure 2, the metabolic model is presented in a simplified way including the cellular compartmentalization, the pathways involved, and the connections between them as well as the exchange reactions with the atmosphere and between day periods.

FIGURE 2

Figure 2 Structure of the diel model presented. Four compartments (chloroplast, mitochondria, peroxisomes, and cytoplasm) and the two phases of the day (light and dark) are considered with exchange reactions including metabolites diffused through leaves (CO₂, O₂, H₂O, light photons) and through the roots (HNO₃ and H₂SO₄). PG represents 2-phospho-glycolate.

The diel model is composed by different organelles including chloroplasts, mitochondria, peroxisomes, and cytoplasm. The proposed model also features plants’ simultaneity of metabolic pathways (i.e., glycolysis in chloroplast and cytoplasm or the folate metabolism in chloroplasts, cytoplasm, and mitochondria) as well as the aforementioned coordination between day phases. A summary of the metabolic model is found in Supplementary S3, and the model in SBML format is accessed in Supplementary S4.

The metabolic network model is mathematically formalized as a constraint-based model, and the fluxes are calculated by applying a flux balance analysis (Cheung et al., 2014):

$\begin{array}{ll}\begin{array}{ll}\begin{array}{c}max\\ v\end{array}\hfill & (v_{biomass}\left[d\right])\hfill \\ subject:\hfill & \begin{array}{l}lb\le Sv\le ub\\ Cv=F\end{array}\hfill \end{array}& (27)\end{array}$

In the FBA formulation in (27), the objective function refers only to the daily biomass production. It is well known for other crop species like Arabidopsis thaliana that biomass growth also takes place during the dark phases of the day (Gomes et al., 2015). However, night metabolism of non-starchy crops like lettuce is still not clear so in the current approach maintenance-associated reactions are limited to night metabolism whereas metabolism associated with light periods of the day concentrates, on top of maintenance, also biosynthetic reactions. Letter υ represents the array of fluxes for each of the reactions of the model. The lower and upper bounds (lb and ub) are fixed only for those fluxes indicated in Table 3. Matrix C contains information regarding the flux ratios specified in Table 3, and matrix F represents the resulting flux. The inequality constraints represented in Table 3 are generated in levels 3 and 2 of the model and are used to feed the FBA. ATP_maintenance and NAD(P)H_maintenance includes the reactions that contribute to the consumption of ATP and the reducing agent for respiration purposes. As suggested in Cheung et al. (2013), this can be achieved including generic ATPase and NAD(P)H oxidase reactions. Finally, the enzyme rates included are those related to constrain the reducing power supply in the cytoplasm. At night, plastidic NADP-malate dehydrogenase and glyceraldehyde 3-phosphate dehydrogenase are downregulated (Miginiac-Maslow and Lancelin, 2002; Mekhalfi et al., 2014). Finally, minimums and maximums for a set of reactions are defined given the presence of thermodynamically infeasible loops when no restrictions are applied. These reactions include the following: PYK, pyruvate kinase; PGM: phosphoglycerate mutase; ENO: enolase; EB1: inorganic pyrophosphatase; EB2: inorganic pyrophosphatase; ACS: acetyl-CoA synthetase; Ser_bio_cl: phosphorylated serine pathway; GOGAT: glutamate synthase; Prot32: 3-mercaptopyruvate sulfurtransferase/cytoplasmic aspartate aminotransferase; OASTL: cysteine synthase; GS: glutamine synthetase. The FBA presented is implemented both in Matlab ® 2021 using the Cobra Toolbox (Becker et al., 2007) and in python 3.0 using Cobrapy (Ebrahim et al., 2013) and can be found in Supplementary S5.

TABLE 3

Table 3 Flux balance analysis equality, inequality, and flow ratio restrictions. [d] and [n] indicate day and night phase period respectively.

2.1.5 Dynamic model

The evolution of the different states of interest through time for either dry biomass, leaf temperature, or gas compositions can be obtained integrating their rates of generation or consumption over time. For dry biomass, this can be done straightforward from Pn, considering that all carbon molecules captured by the plant are fixed into structural biomass:

$\begin{array}{ll}{M}_x=\underset{t=t0}{\overset{t= t_f}{\int }}Pn·BCmol·dt& (28)\end{array}$

For the evolution of leaf temperature, it is necessary to solve an energy balance between the leaf temperature and the environment and to convert energy units to temperature degrees:

$\begin{array}{ll}{T}_l=\underset{t=t0}{\overset{t= t_f}{\int }}\frac{(E_{hs}+E_{hd})A-(E_r+E_{conv}+E_{tr})LA}{C_p\frac{M_x}{DM}}·dt& (29)\end{array}$

Finally, the oxygen concentration can be obtained by solving a mass balance within the growing crop chamber:

$\begin{array}{ll}{O}_2=\underset{t=t0}{\overset{t= t_f}{\int }}\frac{u·(O_2^{in}-O_2)+v_{E{x}_{O_2}}}{\text{V}}·dt& (30)\end{array}$

In (30), gas flow is represented by u, oxygen concentration in the input flow by $O_2^{in}$, and the chamber volume is V.

2.2 Integration of the multilevel mechanistic model with advanced control architectures

Once having defined the modelling strategy for higher crops, the second scope of this study is to integrate the use of metabolic models into an advanced control strategy. In different complex systems ranging from microgrids (Vasquez et al., 2010), life support (Ciurans et al., 2021), or water distribution systems (Ocampo-Martinez et al., 2012) to chemical plants (Scattolini, 2009; Marchetti et al., 2014), the use of advanced control architectures has been proven to be efficient in terms of optimal management and control. Advanced control architectures are characterized by hierarchically distributing management and control functions in different levels. In this study, an adaptation of a common control architecture used in process plants is adapted to control oxygen in a crop-growing chamber. The top layer called steady-state target optimization (SSTO) aims at finding reference values for the controlled and manipulated variables given a specific setpoint through solving a mass balance problem at steady state. The output of SSTO are the controlled and manipulated variables that give the closest estimation of the controlled variables to the setpoint at steady state. This output becomes the input of the following control step in the hierarchy, which is a model-based predictive control (MPC). MPC uses a discretized model of the process to be controlled and aims at finding the sequence of control commands along a prediction horizon that brings the predicted controlled variables the closest to the reference given a set of hard and soft constraints. MPC works based on a rolling-horizon approach, which essentially solves the minimization problem along the defined prediction horizon but sends solely the control command corresponding to the first step. This process is repeated every time the controller is executed (Pannocchia and Bemporad, 2007). In Figure 3, a schematic representation of the advanced control architecture with details on the communication among its levels is presented. In the presented study, the controlled variable is the oxygen concentration in the chamber whereas the manipulated variable is the gas flow.

FIGURE 3

Figure 3 Scheme of the advanced control architecture proposed to integrate metabolic models.

The SSTO solves the following system of equations using the output of FBA, considering oxygen exchange rate (v_ExO2 ) as the generation rate in a steady-state mass balance:

$\begin{array}{ll}{u}_s(O_2^{in}-x_s)+v_{E{x}_{O_2}}=0& (31)\end{array}$

$\begin{array}{ll}{x}_s+d=y_s& (32)\end{array}$

$\begin{array}{ll}{y}_s=y_{ref}& (33)\end{array}$

s.t.

$\begin{array}{ll}{y}_L\le y_s\le y_U& (34)\end{array}$

$\begin{array}{ll}{u}_L\le u_s\le u_U& (35)\end{array}$

As indicated in Figure 3, the output of the SSTO is provided as a reference to the MPC, which in this case is the concentration of oxygen in the chamber (x_s) and the external gas flow (u_s) both at steady state. The internal model at steady state is defined in (3), with the generation term (v_ExO2 ) being the output of the FBA, $O_2^{in}$ the input oxygen concentration, and the internal model prediction defined by x_s. A disturbance (d) is incorporated in (32) to take into account any possible plant-model mismatch or a measured perturbation and must be taken into account for the new reference generation. Hence, the final prediction value (y_s) at steady state considering the presence of any given disturbance will match the process measurement guaranteeing offset-free control (Pannocchia and Bomperad, 2007). The expression to obtain the disturbance is defined in the MPC development hereafter. The technical upper and lower bounds of controlled and manipulated variables are summarised in (34)–(35) and thus do not need to be defined in the MPC. It is then possible to violate constraints on controlled and measured variables during transition states but not at steady state (Marchetti et al., 2014).

The MPC solves a rolling-horizon non-linear optimization problem, taking the output of the SSTO as the reference to track:

$\begin{array}{ll}J=\underset{U}{min }{\lambda }_1\sum _{i=k}^{i=N_p}\frac{{\left|\widehat{x}\left(i|k\right)-x_s\right|}^2}{\Delta x}+{\lambda }_2\sum _{i=k}^{i=N_c}\frac{{\left|u(i|k)-u_s\right|}^2}{\Delta u}+{\lambda }_3\sum _{i=k}^{i=N_c}\frac{{\left|u(i|k)-u(i-1|k)\right|}^2}{\Delta u}& (36)\end{array}$

s.t.

$\begin{array}{ll}\widehat{x}(i+1|k)=\widehat{x}(i|k)+U(i|k)\left(O_2^{in}-\widehat{x}(i|k)\right)\frac{T_s}{V}+v_{E{x}_{O_2}}\frac{T_s}{V}& (37)\end{array}$

$\begin{array}{ll}\widehat{y}(i|k)=\widehat{x}(i|k)+\text{d(}i|k)& (38)\end{array}$

$\begin{array}{ll}\text{d(}i|k)=y_p-\widehat{x}(2|k-1)& (39)\end{array}$

$\begin{array}{ll}\text{d(}i+1|k)=\text{d(}i|k)& (40)\end{array}$

$\begin{array}{ll}U=\left[\begin{array}{c}U(i\mid k)\\ ⋮\\ \text{U}(N_c\mid \text{k})\end{array}\right]& (41)\end{array}$

The cost function in (36) includes penalization terms to the deviation of the internal prediction ($\widehat{x}$) and the control command (u) from the concentration (x_s) and the gas flow (u_s) references generated in the SSTO. A third penalization term to the rate of change of the manipulated variable is also included in (36) aimed to adjust the speed of the controller response. All three penalization terms are normalized using the range of possible maximum and minimum values for both controlled and manipulated variables (defined by Δx and Δu) and are subject to a scaling factor (λ). Prediction and control horizons are represented by N_p and N_c, respectively. Constraints in (37) reflect the dynamics of the growing crop chamber with T_s as the sample time of the MPC. The internal model is initialized with the current process measurement (y_p). In (39), the disturbance is integrated to the internal model prediction as similarly done in the SSTO in (32). This disturbance is estimated at each sampling time and is defined as the difference between the process measurement (y_p) and the first step prediction of the previous MPC execution ($\widehat{y}\left(2|k-1\right)$) as stated in Tatjewski (2017). It is assumed in (40) the disturbance to be constant through the whole prediction horizon. The output matrix of gas flows (U) is expressed in (41), and only the control action for the first step of the control horizon is sent to the control actuators until the next SSTO and MPC execution. This control strategy is implemented in Matlab® 2021 using the Optimization Toolbox for solving non-linear programming problems and Cobra Toolbox 2022 for the FBA resolution (Becker et al., 2007).

2.3 Simulation scenarios

Two simulation packages are presented: first, the results of the multilevel model presented in Section 2.1 and their validation with experimental data, and second, a dynamic simulation presenting the response of the control architecture presented in Section 2.2 under different perturbations.

2.3.1 Simulation conditions for multilevel model validation

CO₂-response curves were generated using the fixed light intensity indicated in Table 4 and the following range of internal CO₂ values in μmol CO₂/mol: 100, 225, 300, 450, 600, 850, 1,000, 1,100. Light-response curves were generated using the fixed internal carbon dioxide concentration indicated in Table 4 and the following range of light intensity values in μmol/m²/s: 100, 200, 350, 500, 600, 800, 1,000, 1,200. Considering the conditions listed in Table 4 and the model Equations (1)–(26), it is possible to retrieve the inequality constraints described in Table 3 and thus the FBA can be resolved. Results are represented in Section 3.1. The distribution of fluxes is analysed for both light and dark metabolism using atmospheric conditions for CO₂ which are 400 ppm and a light intensity of 400 µmole/m²/s. Results are graphically represented in Section 3.2.

TABLE 4

Table 4 List of operating parameters used in FBA simulation.

2.3.2 Advanced control architecture configuration

In Table 5, the controller specifications and parameter values are indicated. On top of the control objective defined in Section 2.2, atmospheric CO₂ is controlled at 800 ppm with external addition of pure CO₂.

TABLE 5

Table 5 Controller specifications including SSTO and MPC algorithms.

Worthy of note is that the scaling factors are used to promote the control of the system close to the setpoint but at the expense of having a more aggressive control. Prediction and control horizons are important tuning parameters of the controller increasing its sensitivity but also the computational cost of the calculation. Finally, the sampling time also affects the control performance. For slow systems like a plant cultivar, sample times should not be too short because the prediction would not have enough perspective to take correct decisions. Tuning model-based predictive controllers is critical to achieving a good process operation and represents a trade-off between the expected performance and the controller and system capabilities. The metabolic-based control architecture represented in Figure 3 is tested in a 24-h dynamic simulation using the operating conditions listed in Table 4.

3 Results

This section is divided into two parts: the first dedicated to present the output of the multilevel model and its validation with experimental data and the second related to the integration of that modelling approach containing metabolic information into an advanced control architecture.

3.1 Validating photosynthesis rates

The results of the model introduced in Section 2.1 are represented in Figure 4, where modelled and experimental results for mature leaves (28 days after transplanting) are compared.

FIGURE 4

Figure 4 Model validation through the comparison of the photosynthesis rate expressed as net carbon assimilation (Pn) as a function of internal carbon concentration and light irradiance. In (A), a fixed light intensity of 800 μmol/m²/s is set; in (B), a fixed internal carbon concentration of 400 μmol/m²/s is set. Experimental data are based on Zhou et al. (2020).

Figure 4A represents the RuBisCo saturation curve showing a fast rate of change in the smaller range of internal carbon dioxide concentrations, a pattern which is reproduced in Figure 4B which shows the light electron transport chain saturation curve. The maximum net photosynthesis rate achieved in the light-response curve is lower than the maximum achieved in the CO₂-response curve because in the former, the internal CO₂ concentration used for the simulation is 400 μmol CO₂/mol reaching the expected photosynthesis rates if compared to Figure 4A. Finally, the output of the model is comparable to the reported experimental results by Zhou et al. (2020) under the same operating conditions. Overall, the error observed in the modelled results in relation to the experimental values is higher at low internal carbon concentrations and at low light intensities.

3.2 Distribution of metabolic fluxes using a metabolic diel model

In this section, the flux distribution of day and night metabolites obtained after the resolution of the FBA introduced in (27) are presented. The parameters used for the simulation are those provided in Table 3.

3.2.1 Day flux distribution

The day flux distribution is presented in Figure 5. The central carbon metabolism of plants in light conditions is well represented in this model with the main fluxes located in the Calvin cycle-associated reactions. The results indicate that, as has been extensively studied and published (Michelet et al., 2013; Tan and Cheung, 2020), the flow through the Calvin cycle generates triose phosphate (g3p) from ribulose 1,5-biphosphate (RuBP), consuming part of the reducing power and ATP molecules synthesized during the light electron transport chain. Triose phosphate is used to feed the rest of the Calvin cycle machinery aimed at restoring the ribulose 1,5-biphosphate while it is also partially used to generate photosynthetic end products (McClain and Sharkey, 2019).

FIGURE 5

Figure 5 Flow distribution of the central carbon metabolism during L. sativa grown during the light photoperiod. List of abbreviated enzymes or enzyme reactions. In chloroplast: F_LETC, flow of light electron transport chain; G6Pase, glucose-6-phosphatase; PGI, glucose-6-phosphate isomerase; MDH, malate dehydrogenase; FBPase, fructose-1,6-biphosphatase; PFK, phosphofruktokinase-1; ALD, aldolase; TKT, transketolase; TKT2, transketolase 2; FBA, fructose-1,6-biphosphate aldolase; TPI, triose-phosphate isomerase; RPI, ribose-5-phosphate isomerase; RPE, ribulose-phosphate 3-epimerase; GAPDHy, NADP-glyceraldehyde-3-phosphate dehydrogenase; PGK, phosphoglycerate kinase; PRK, phosphoribulokinase; V_c, RuBisCo carboxylation; V_o, RuBisCo oxgygenation; Glyk, D-glycerate 3-kinase; PGP, phosphoglycolate phosphatase. In peroxisome: GOXp, glycolate oxidase; SGTp, serine-glyoxylate transaminase; GTp, glycine transaminase; HPRp, hydroxypyruvate reductase; MDHp, malate dehydrogenase. In cytoplasm: PGMc, phosphoglycerate mutase; ENOc, enolase; PYKc, pyruvate kinase; PGKc, phosphoglycerate kinase; GAPDHc, NAD-glycerate-3-phosphate dehydrogenase; MDHc, NADH malate dehydrogenase; ICDHyc, NADP-based isocitrate dehydrogenase; ALDc, aldolase; PFK, phosphofructokinase; TALc, transaldolase. In mitochondria: ACONTm, aconitase; MDHm, malate dehydrogenase; GCSm, glycine cleavage system; SHMTm, serine hydroxymethyltransferase; CSm, citrate synthase; PCm, pyruvate carboxylase; MP2, mitochondrial phosphorylation 2. Carbohydrate reactions represent a set of lumped reactions related to carbohydrate metabolism.

As previously stated, higher plants store carbon during the light phase of the day, to be used during night respiration and fuel maintenance processes. For non-starchy crops, even though they can still generate starch, most of the carbon fixed during the light phase is stored as soluble sugars or organic acids. In this study, the sugar molecules stored and mobilized between light and dark periods of the day have not been restricted and, as Figure 5 suggests, sucrose, fructose, citrate, malate, and fumarate are the metabolites used for carbon exchange. This modelled result fits well with the reported experimental concentrations of sugars in lettuce at harvest, with glucose, sucrose, and fructose being the main carbohydrates found for carbon exchange between phases of the day (Chen et al., 2019) and also predicting with accuracy the role of malate accumulation during the light phase of the day in vacuoles for its use in the dark (Lee et al., 2021). Not all carbon compounds mobilized from light to night metabolism are consumed during the latter. Therefore, some carbon intermediates need to be exported from night to light metabolism too. Specifically, citrate imported from dark periods is used in the light phase of the day to generate oxoglutarate which is important for the nitrogen assimilation mechanism and for the synthesis of nitrogen-rich amino acids. The citrate cycle and its interactions with amino acid biosynthesis are well covered by the presented model, both suffering a flux reduction when nitrate uptake is limited (Morcuende et al., 1998; Weiwei Zhou et al., 2021).

One of the critical phenomena of plant photosynthetic cell metabolism is the coordination of photosynthesis and respiration, which essentially determines how and where energy carrier molecules (ATP and NAD(P)H) are produced. Most of the ATP and NAD(P)H used for catabolic reactions are produced in the light electron transport chain in chloroplasts for amino acid and lipid production. Part of the ATP synthesized in chloroplasts is exported to the cytosol through the 3PG-G3P shuttle, satisfying the ATP demand together with ATP exported from mitochondria (Gakière et al., 2018; Shameer et al., 2019). The resulting metabolic network shows the mechanisms of redox power balancing in the different organelles of photosynthetic leaves enabled by the metabolite shuttles represented by the malate/oxalacetate and the triose phosphate/3-phosphoglycerate shuttle in Figure 5 (and the glutamate/2-oxoglutarate and the malate/aspartate shuttle, not represented) (Taniguchi and Miyake, 2012). Around half the NAD(P)H generated in mitochondria comes from the glycine decarboxylation, which in turn generates the serine used in the serine-glyoxylate aminotransferase (SGT) in the peroxisome. It has been reported that malate dehydrogenase seems to regulate the reducing power in mitochondria based on the reduction state of the cells, removing and restoring NAD(P)H at low- and high-light conditions responding to changes in the photosynthesis rate (Bykova et al., 2014; Schertl and Braun, 2014). This is validated in the presented fluxome, where mitochondrial malate dehydrogenase (MDHm) removes excess NAD(P)H produced through the glycine cleavage system (GCS, also known as glycine decarboxylase system), the main contributor of redox power in the mitochondria. The metabolic flux distribution represented in Figure 5 also shows that the TCA cycle is not complete during light photoperiods. This is mainly because most TCA intermediates are dedicated to anabolic reactions and pyruvate dehydrogenase (PDH) is photo-inhibited (Schertl and Braun, 2014).

3.2.2 Night flux distribution

During the dark phase of the day in Figure 6, ATP can only be produced in mitochondria with NAD(P)H being the electron donor and thus completely modifying the flow distribution within the cell. Autotrophic organisms like plants use the sugars generated and stored during the day to feed TCA, which is cyclic in dark conditions. In lettuce and other non-starchy vegetables, instead of mobilizing starch, glucose, sucrose, or fructose is broken down to pyruvate at very similar proportions (Chen et al., 2019) even though in this study the proportion of soluble sugar utilization has not been constrained. The pyruvate produced in the glycolysis is then transferred to the mitochondria to regenerate the reducing power needed to fuel the mitochondrial respiration. In dark conditions, the exchange rates are completely opposite to those observed in light conditions, with carbon dioxide and water being released and oxygen consumed. In a way, respiration and photosynthesis are opposite processes, but complementary as demonstrated during crop growth. The governing rule of night metabolism is represented in Table 3 in the ratio of consumed CO₂ during the day phase over produced CO₂ due to night respiration (Ex_CO2[d]:Ex_CO2[n]) which sets the night carbon conversion.

FIGURE 6

Figure 6 Flow distribution of the central carbon metabolism during L. sativa dark phase. List of abbreviated enzymes or enzyme reactions. In cytoplasm: MDHc, malate dehydrogenase; ICDHyc, NADPH isocitrate dehydrogenase; ACONT, aconitase; MEc, malic enzyme; CP1, glucose isomerase; CP2, sucrose-6-phosphate synthase. In mitochondria: ACONTm, aconitase; PDHm, pyruvate dehydrogenase; PCm, pyruvate carboxylase; CSm, citrate synthase; ICDHm, isocitrate dehydrogenase; AKGDm, α-oxoglutarate dehydrogenase; SUCOASm, succinyl-CoA synthetase; SUCD1m, succinate dehydrogenase complex; MDHm, malate dehydrogenase; MP1/2, mitochondrial phosphorylation 1/2.

3.2.3 Sensitivity analysis

FBA is a powerful tool but strongly affected by its mathematical formalization and particularly by the flux boundary choice (Raposo et al., 2020; Nobile et al., 2021). A sensitivity analysis on the selected boundaries is very useful for detecting those constraints (either hard or soft constraints) that generate the highest impact on the flux distribution and thus requires a thorough parameter identification. However, sensitivity analysis can also be used to detect those fluxes with the highest variability in relation to a given boundary, metabolic model shortcomings, or fluxes that are invariable to the boundaries. A local sensitivity analysis is presented, with a focus on the flux ratios defined in Table 3, which are either defined empirically or based on literature and may be prone to uncertainty. The sensitivity analysis has also been explored for the light irradiance. Details on the ranges of flux ratios explored for this analysis are presented in Table 6. The sensitivity analysis has been based on the variability of the fluxes along the range of flux ratios allowed. Such variability is expressed as the normalized summation of the slopes in each of the steps of the range of flux ratios explored:

$\begin{array}{ll}F{V}_{j,z}= \sum _{i=1}^{i=6}\frac{\frac{{\nu }_{i+1,j,z}-{\nu }_{i,j,z}}{\text{Δ}{r}_z}}{\text{Δ}{v}_{j,z}}& (42)\end{array}$

$\begin{array}{ll}\text{Δ}{v}_{j,z}= \max ({\nu }_{j,z})-\min ({\nu }_{j,z})& (43)\end{array}$

TABLE 6

Table 6 Range of values explored in sensitivity analysis.

In (42), the i index refers to the discrete steps in which the range of flux ratios has been divided (a total of six steps for each flux ratio) whereas index j and z refer to the reaction and the flux ratio, respectively. The parameter Δr_z is the increment in the flux ratio expressed as a fraction and equivalent to 0.1 for all cases. The denominator Δv_j,z is used to normalize the summation of slopes, and it represents the range of values for a given reaction j. The sensitivity analysis has been carried for light and night metabolism and is represented in Figure 7 considering an atmospheric concentration of CO₂ of 1,000 ppm in order not to limit by carbon substrate and a nominal light intensity of 400 µmol/m²/s.

FIGURE 7

Figure 7 Normalized sensitivity analysis result of a selection of day (A) and night (B) metabolic fluxes (B) according to a given range of ratios and constraints explored.

Considering the light metabolism sensitivity analysis (Figure 7A), the light phase of day-time mitochondrial respiration is never active. In the predicted fluxome for all ranges of flux ratios applied, the ATP production comes mainly from chloroplast activity. The main mitochondrial NADH producer in light conditions is the GDC system, being consumed mainly by the malate dehydrogenase activity and not by the electron transport chain (Bykova et al., 2014). The lack of the electron transport chain in day metabolism is due to the lack of maintenance reactions during the light period. When light is increased, the reactions associated with the regeneration of RuBP are activated. This is the case of phosphoribulokinase (PRK), which restores the ribulose 1,5-biphosphate. An increase in the triose phosphate export activity (G3P_trans) is also detected responding to the increase in the carbon fixation. The mitochondrial malate dehydrogenase (MDHm) is also augmented in the first place to respond to the increased anabolic demands and in the second place given the increased flux through photorespiration that triggers the GCS and the consequent increased consumption of NADH through MDHm to keep the redox balance in mitochondria (Schertl and Braun, 2014). The sensitivity analysis of the day metabolism also highlights that, when the night respiration flux in relation to the daily carbon fixation is increased relative to the nominal value of 0.25, the associated photosynthetic reactions (PRK and G3P_trans) are reduced.

Regarding night metabolism (Figure 7B), the variability of the represented reactions is higher than in day metabolism. When light intensity (I_u) is increased, night mitochondrial activity (MP1_n and MP2_n) is also increased to respond to the maintenance reaction demands (Frantz and Bugbee, 2005). The activity of glyceraldehyde 3-phosphate dehydrogenase (GAPDc_n) responds to the increased flux towards glycolytic pathways to process the carbon compounds converted during night metabolic activities (Schneider et al., 2018; Gaude et al., 2018). In the dark phase, with the cyclic TCA cycle re-stored, the excess TCA intermediates and soluble sugars not used to fuel dark metabolism are stored and reused in light metabolism for nitrogen assimilation and biosynthesis reactions (Popova and Pinheiro De Carvalho, 1998; Igamberdiev and Eprintsev, 2016). In this way, the pool of organic carbon compounds and reducing equivalents is managed by the plants as a response to variations in the day carbon fixation efficiency and night respiration activities. GAPDc_n is also positively affected when the ratio of respiratory produced over photosynthetically consumed CO₂ is increased (Ex_CO2[n]:Ex_CO2[d]) with the availability of citrate and the usage of other TCA intermediates in light photoperiods reduced given the night metabolic increased activities as suggested by the citrate import flux. Therefore, an increased night CO₂ release is corresponded by an increased glycolytic activity represented by GAPDC_n, also impacting the TCA intermediates and soluble sugar accumulation at night.

3.3 Testing the integration of multilevel model approach to advanced control architectures in dynamic simulations

The control strategy presented based on a new approach combining the prediction capacity of MPC and the integration of constraint-based metabolic modelling has been tested in a 24-h dynamic simulation. To do that, different perturbations have been included at different points of the simulation to test the resiliency of the proposed control approach its capacity to overcome common control challenges and to analyse its versatility:

● At 5 h of simulation time, a plant-model mismatch has been introduced by adding a perturbation in the form of a multiplication factor to the gas flows in (31):

$\begin{array}{ll}2u_s(O_2^{in}-x_s)+v_{E{x}_{O_2}}=0& (44)\end{array}$

$\begin{array}{ll}\widehat{x}(i+1|k)=\widehat{x}(i|k)+2U(i|k)\left(O_2^{in}-\widehat{x}(i|k)\right)\frac{T_s}{V}+v_{E{x}_{O_2}}\frac{T_s}{V}& (45)\end{array}$

● • At 10 h of simulation time, a perturbation in the process output has been introduced by adding a sudden decrease in the oxygen concentration from the measurement to 20.8%.

● • At 16 h of simulation time, a change in the oxygen setpoint from 21% to 21.2% has been included.

The results of this simulation schedule are presented in Figure 8. It is observed that when a plant-model mismatch is deliberately included, the controller can keep minimising the offset of the O₂ measurement in relation to the reference (Figure 8A). Similarly, when an abrupt perturbation is added at 10 h of simulation, the gas flow is stopped to restore the oxygen concentration rapidly. Finally, when the reference is modified from 21% to 21.2%, the gas flow is also reduced in order to accumulate oxygen in the growing chamber and reduce the tracking error (Figure 8B). The control of CO₂ is achieved by externally injection of pure CO₂ (Figure 8C).

FIGURE 8

Figure 8 Dynamic control performance under introducing different perturbations at 5, 10, and 16 hours as described in the text. (A) Tracking of the Oxygen concentration; (B) Evolution of the gas flow as a manipulated variable; (C) Dynamics of concentration and external addition of Carbon dioxide.

The controller is thus demonstrated to be resilient and smooth to overcome any of the perturbations applied as well as on the nominal operation. This performance is due to the reliable metabolic-based model introduced in the MPC and largely due to the integration of a disturbance specially to achieve offset-free control. In Figure 9, the comparison of the controller performance with and without disturbance integration is represented. It can be observed that under different perturbations, when disturbance is integrated in the internal model of both the SSTO and the MPC, offset is reduced and thus the reference tracking is improved. Nevertheless, not all disturbances can be rejected using a constant state disturbance prediction as defined in (32), especially in scenarios where plant-model mismatches are bigger, that is, when models are less reliable than the one presented in this study. In the case of higher model disagreement, disturbance integration approaches should be considered like integrating the error in the multistep model prediction (Tian et al., 2007) or integrating a moving horizon estimation (MHE) and MPC to estimate uncertainty parameters and to include them in the MPC algorithm (Huang et al., 2010). Disturbance rejection especially in the situation of plant-model mismatch has been highly analysed within MPC development, and certainly, the future metabolic-based controllers will need to deal with a wide diversity of model typologies from simplified and surrogated models to complex genome-scale metabolic models.

FIGURE 9

Figure 9 Comparison of the controller performance in terms of reference tracking between controllers including or not disturbance rejection in the internal SSTO and MPC models.

4 Conclusions

Modelling higher plants has major challenges to deal with. These challenges include huge metabolic changes associated with the light and night photoperiods, substrate partitioning given the heterogenic requirements of the different tissues present in higher plants, organelle coordination, complex morphologies that condition the interaction with the environment, and many other phenomena still not fully understood. In this study, a multilevel model has been designed with the main mechanistic phenomena that drive crop growth distributed into different levels in decreasing order of scale length. The output of the mechanistic multilevel model has been connected to a constraint-based metabolic model providing information of high interest about cell metabolism. The presented multilevel model offers the advantage to merge all the available information related to plant growth in a structured way ensuring the solution found to be feasible for all the phenomena described in the different layers in the model hierarchy, from light reception and biochemical conversion down to specific metabolic pathways, a feasibility that cannot be granted in only mechanistic or metabolic-based models using stand-alone modelling strategies. This method has been validated with experimental data, integrated, and tested in a novel advanced control strategy with promising results. Computational capabilities are no longer a major constraint, making it possible to contemplate the design of control strategies that can integrate more information about the system under operation than currently done. In this study, the focus has been placed in the use of metabolic information embedded in a model-predictive control providing promising results in a dynamic simulation of a growing crop chamber with a range of applications going from agriculture to life support systems. It is in the framework of the MELiSSA project and by extension to the field of fully regenerative life support systems where structural multilevel models, integrating from physical to metabolic information, emerges as an opportunity to design model-based predictive controller techniques. These advanced control strategies should be further explored with different levels of metabolic complexity, in different control formulations, and applied to different biobased processes as they may contribute to improving the overall performance through exploiting the increasingly available metabolic information.

Data availability statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

Author contributions

CC is the main writer of the manuscript as well as the one responsible for script development, simulation, and main analysis and interpretation of the results. IM and IMdM are the main contributors to the constraint-based metabolic modelling. IMdM contributed to the design, validation, and interpretation of the results. JM is the main contributor and is responsible for the design of a control architecture adapted to metabolic networks. CG and FG have contributed to the design of the mechanistic multilevel model, specially concerning the physical and enzymatic equations, and they have supervised the overall text. CG is the responsible of the original lettuce metabolic model development. All authors contributed to the article and approved the submitted version.

Funding

MELiSSA is an international consortium of 15 partners led by the European Space Agency. Its activities are governed by a Memorandum of Understanding (ESA 4000100293/10/NL/PA). The MELiSSA Pilot Plant is funded from ESA contributions from Spain (main contributor), Belgium, France, Italy, and Norway, under Frame Contract C4000109802/13/NL/CP. Cofunding from Ministerio de Ciencia, Innovación y Universidades, Generalitat de Catalunya, and Universitat Autònoma de Barcelona is also acknowledged. Carles Ciurans is a PhD fellow from the POMP Program (MELiSSA Foundation) of ESA. This work was supported by VILLUM FONDEN under the VILLUM Investigator Grant (no. 25920): Center for Research on Microgrids (CROM). Igor Marin de Mas and Ivan Martínez-Mongue would like to thank The Novo Nordisk Foundation (NNF Grant numbers: NNF10CC1016517 and NNF14OC0009473).

Acknowledgments

The authors wish to thank the partners of the MELiSSA Consortium for carefully reviewing this publication.

Conflict of interest

Author IMdM and IM are employed by Novo Nordisk Foundation Center for Sustainability.

The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fpls.2022.970410/full#supplementary-material

References

Amitrano, C., Chirico, G. B., Pascale, S., Rouphael, Y., Micco, V. (2020). Crop management in controlled environment agriculture (CEA) systems using predictive mathematical models. Sens. (Switzerland) 20 (11), 1-19. doi: 10.3390/s20113110

CrossRef Full Text | Google Scholar

Arnold, A., Nikoloski, Z. (2011). A quantitative comparison of Calvin-Benson cycle models. Trends Plant Sci. 16 (12), 676–683. doi: 10.1016/j.tplants.2011.09.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Becker, S. A., Feist, A. M., Mo, M. L., Hannum, G., Palsson, B., Herrgard, M. J. (2007). Quantitative prediction of cellular metabolism with constraint-based models: The COBRA toolbox. Nat. Protoc. 2 (3), 727–738. doi: 10.1038/nprot.2007.99

PubMed Abstract | CrossRef Full Text | Google Scholar

Boscheri, G., Kacira, M., Patterson, L., Giacomelli, G., Sadler, P., Furfaro, R. (2012). Modified energy cascade model adapted for a multicrop lunar greenhouse prototype. Adv. Space Res. 50 (7), 941–951. doi: 10.1016/j.asr.2012.05.025

CrossRef Full Text | Google Scholar

Bykova, N. V., Møller, I. M., Gardeström, P., Igamberdiev, A. U. (2014). The function of glycine decarboxylase complex is optimized to maintain high photorespiratory flux via buffering of its reaction products. Mitochondrion 19, 357–364. doi: 10.1016/j.mito.2014.01.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, X., Wang, L.c., Li, T., Yang, Q.c., Guo, W.z. (2019). Sugar accumulation and growth of lettuce exposed to different lighting modes of red and blue LED light. Sci. Rep. 9 (1), 1–10. doi: 10.1038/s41598-019-43498-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheung, C. Y. M., Williams, T. C. R., Poolman, M. G., Fell, D. A., Ratcliffe, R. G., Sweetlove, L. J. (2013). A method for accounting for maintenance costs in flux balance analysis improves the prediction of plant cell metabolic phenotypes under stress conditions. Plant J. 75 (6), 1050–1061. doi: 10.1111/tpj.12252

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheung, C. Y M., Poolman, M. G., Fell, D. A., George Ratcliffe, R., Sweetlove, L. J. (2014). A diel flux balance model captures interactions between light and dark metabolism during day-night cycles in C3 and crassulacean acid metabolism leaves. Plant Physiol. 165 (2), 917–929. doi: 10.1104/pp.113.234468

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciurans, C., Bazmohammadi, N., Vasquez, J. C., Dussap, G., Guerrero, J. M., Godia, F. (2021). Hierarchical control of space closed ecosystems – expanding microgrid concepts to bioastronautics. IEEE Ind. Electron. Mag. 15 (2), 2–13. doi: 10.1109/MIE.2020.3026828

CrossRef Full Text | Google Scholar

Corpas, F. J., Barroso, J. B. (2014). NADPH-generating dehydrogenases: Their role in the mechanism of protection against nitro-oxidative stress induced by adverse environmental conditions. Front. Environ. Sci. 2. doi: 10.3389/fenvs.2014.00055

PubMed Abstract | CrossRef Full Text | Google Scholar

Dong, C., Fu, Y., Xie, B., Wang, M., Liu, H. (2017). Element cycling and energy flux responses in ecosystem simulations conducted at the Chinese lunar palace-1. Astrobiology 17 (1), 78–86. doi: 10.1089/ast.2016.1466

PubMed Abstract | CrossRef Full Text | Google Scholar

Ebrahim, A., Lerman, J. A., Palsson, B. O., Hyduke, D. R. (2013). COBRApy: COnstraints-based reconstruction and analysis for Python. BMC Syst. Biol. 7, 1-6. doi: 10.1186/1752-0509-7-74

PubMed Abstract | CrossRef Full Text | Google Scholar

Eckart, P. (1995). Spaceflight life support and biospherics. (Dordrecht: Springer), 1–444.

Google Scholar

Farquhar, G. D., Caemmerer, S.v., Berry, J. A. (1980). A biochemical model of photosynthetic CO₂ assimilation in leaves of C₃ species. Planta 149 (1), 78–90. doi: 10.1007/BF00386231

PubMed Abstract | CrossRef Full Text | Google Scholar

Frantz, J. M., Bugbee, B. (2005). Acclimation of plant populations to shade: Photosynthesis, respiration, and carbon use efficiency. J. Am. Soc. Hortic. Sci. 130 (6), 918–927. doi: 10.21273/jashs.130.6.918

CrossRef Full Text | Google Scholar

Fulget, N., Poughon, L., Richalet, J., Lasseur, C. (1999). MELISSA : Global control strategy of the artificial ecosystem by using first principles models of the compartments. Adv. Space Res. 24 (3), 397–405. doi: 10.1016/S0273-1177(99)00490-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Gakière, B., Hao, J., Bont, L., Pétriacq, P., Nunes-Nesi, A., Fernie, A. R. (2018). NAD + biosynthesis and signaling in plants. Crit. Rev. Plant Sci. 37 (4), 259–307. doi: 10.1080/07352689.2018.1505591

CrossRef Full Text | Google Scholar

Gaude, E., Schmidt, C., Gammage, P. A., Dugourd, A., Blacker, T., Chew, S. P., et al. (2018). NADH shuttling couples cytosolic reductive carboxylation of glutamine with glycolysis in cells with mitochondrial dysfunction. Mol. Cell 69 (4), 581–593.e7. doi: 10.1016/j.molcel.2018.01.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Gitelson, J., Lisovsky, G. (2002). Man-made closed ecological systems (London: CRC Press).

Google Scholar

Gomes, C., Dal’molin, O., Quek, L.-E., Saa, P. A., Nielsen, L. K., Nikoloski, Z., et al. (2015). A multi-tissue genome-scale metabolic modeling framework for the analysis of whole plant systems. Front. Plant Sci 6. doi: 10.3389/fpls.2015.00004

CrossRef Full Text | Google Scholar

Harley, P. C., Thomas, R. B., Reynolads, J. F., Strain, B. R. (1992). Modelling photosynthesis of cotton grown in elevated CO2. Plant Cell Environ. 15 (3), 271–282. doi: 10.1111/j.1365-3040.1992.tb00974.x

CrossRef Full Text | Google Scholar

Huang, R., Biegler, L. T., Patwardhan, S. C. (2010). Fast offset-free nonlinear model predictive control based on moving horizon estimation. Ind. Eng. Chem. Res. 49 (17), 7882–7890. doi: 10.1021/ie901945y

CrossRef Full Text | Google Scholar

Igamberdiev, A. U., Eprintsev, A. T. (2016). Organic acids: The pools of fixed carbon involved in redox regulation and energy balance in higher plants. Front. Plant Sci. 7. doi: 10.3389/fpls.2016.01042

CrossRef Full Text | Google Scholar

Lee, C. P., Elsässer, M., Fuchs, P., Fenske, R., Schwarzländer, M., Millar, A. H. (2021). The versatility of plant organic acid metabolism in leaves is underpinned by mitochondrial malate-citrate exchange. Plant Cell 33 (12), 3700–3720. doi: 10.1093/plcell/koab223

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., van Iersel, M. W. (2021). Photosynthetic physiology of blue, green, and red light: Light intensity effects and underlying mechanisms. Front. Plant Sci. 12. doi: 10.3389/fpls.2021.619987

CrossRef Full Text | Google Scholar

Louarn, G., Song, Y. (2020). Two decades of functional-structural plant modelling: Now addressing fundamental questions in systems biology and predictive ecology. Ann. Bot. 126 (4), 501–509. doi: 10.1093/aob/mcaa143

PubMed Abstract | CrossRef Full Text | Google Scholar

Marchetti, A. G., Ferramosca, A., González, A. H. (2014). Steady-state target optimization designs for integrating real-time optimization and model predictive control. J. Process Control 24 (1), 129–145. doi: 10.1016/j.jprocont.2013.11.004

CrossRef Full Text | Google Scholar

Marshall-Colon, A., Long, S. P., Allen, D. K., Allen, G., Beard, D. A., Benes, B., et al. (2017). Crops in silico: Generating virtual crops using an integrative and multi-scale modeling platform. Front. Plant Sci. 8. doi: 10.3389/fpls.2017.00786

CrossRef Full Text | Google Scholar

McClain, A. M., Sharkey, T. D. (2019). Triose phosphate utilization and beyond: From photosynthesis to end product synthesis. J. Exp. Bot. 70 (6), 1755–1766. doi: 10.1093/jxb/erz058

PubMed Abstract | CrossRef Full Text | Google Scholar

Mekhalfi, M., Puppo, C., Avilan, L., Lebrun, R., Mansuelle, P., Maberly, S. C., et al. (2014). Glyceraldehyde-3-Phosphate dehydrogenase is regulated by ferredoxin-NADP reductase in the diatom asterionella Formosa. New Phytol. 203 (2), 414–423. doi: 10.1111/nph.12820

PubMed Abstract | CrossRef Full Text | Google Scholar

Michelet, L., Zaffagnini, M., Morisse, S., Sparla, F., Pérez-Pérez, M. E., Francia, F., et al. (2013). Redox regulation of the Calvin-Benson cycle: Something old, something new. Front. Plant Sci. 4. doi: 10.3389/fpls.2013.00470

PubMed Abstract | CrossRef Full Text | Google Scholar

Miginiac-Maslow, M., Lancelin, J. M. (2002). Intrasteric inhibition in redox signalling: Light activation of NADP-malate dehydrogenase. Photosyn. Res 72 (1), 1-12. doi: 10.1023/A:1016099228450

CrossRef Full Text | Google Scholar

Morcuende, R., Krapp, A., Hurry, V., Stitt, M. (1998). Sucrose-feeding leads to increased rates of nitrate assimilation, increased rates of a-oxoglutarate synthesis, and increased synthesis of a wide spectrum of amino acids in tobacco leaves. Planta 206, 394–409. doi: 10.1007/s004250050415

CrossRef Full Text | Google Scholar

Morgan, J. A., Rhodes, D. (2002). Mathematical modeling of plant metabolic pathways. Metab. Eng. 4 (1), 80–89. doi: 10.1006/mben.2001.0211

PubMed Abstract | CrossRef Full Text | Google Scholar

Nikolov, N. T., Massman, W. J., Schoettle, A. W. (1995). Coupling biochemical and biophysical processes at the leaf level: An equilibrium photosynthesis model for leaves of C3 plants. Ecol. Model. 80 (2–3), 205–235. doi: 10.1016/0304-3800(94)00072-P

CrossRef Full Text | Google Scholar

Nobile, M. S., Coelho, V., Pescini, D., Damiani, C. (2021). Accelerated global sensitivity analysis of genome-wide constraint-based metabolic models. BMC Bioinf. 22, 1-17. doi: 10.1186/s12859-021-04002-0

CrossRef Full Text | Google Scholar

Ocampo-Martinez, C., Barcelli, D., Puig, V., Bemporad, A. (2012). Hierarchical and decentralised model predictive control of drinking water networks: Application to Barcelona case study. IET Control Theory Appl. 6 (1), 62–71. doi: 10.1049/iet-cta.2010.0737

CrossRef Full Text | Google Scholar

Pannocchia, G., Bemporad, A. (2007). Combined design of disturbance model and observer for offset-free model predictive control. IEEE Trans. Automat. Control 52 (5), 1109–1113. doi: 10.1109/TAC.2007.899096

CrossRef Full Text | Google Scholar

Peng, B., Guan, K., Tang, J., Ainsworth, E. A., Asseng, S., Bernacchi, C. J., et al. (2020). Towards a multiscale crop modelling framework for climate change adaptation assessment. Nat. Plants 6 (4), 338–348. doi: 10.1038/s41477-020-0625-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Popova, T. N., Pinheiro De Carvalho, M. A. A. (1998). Citrate and isocitrate in plant metabolism. Biochim. Biophys. Acta 1364, 307–325. doi: 10.1016/S0005-2728(98)00008-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Poughon, L., Farges, B., Dussap, C. G., Godia, F., Lasseur, C. (2009). Simulation of the MELiSSA closed loop system as a tool to define its integration strategy. Adv. Space Res. 44 (12), 1392–1403. doi: 10.1016/j.asr.2009.07.021

CrossRef Full Text | Google Scholar

Poulet, L. (2018). Developpement de modéles physiques pour comprendre la croissance des plantes en environment de gravité réduite pour des applications dans Les systémes suuport-vie (Clermont-Ferrand: Université Clermont-Auvergne).

Google Scholar

Poulet, L., Dussap, C. G., Fontaine, J. P. (2020). Development of a mechanistic model of leaf surface gas exchange coupling mass and energy balances for life-support systems applications. Acta Astronaut. 175, 517–530. doi: 10.1016/j.actaastro.2020.03.048

CrossRef Full Text | Google Scholar

Poulet, L., Fontaine, J. P., Dussap, C. G. (2018). A physical modeling approach for higher plant growth in reduced gravity environments. Astrobiology 18 (9), 1093–1100. doi: 10.1089/ast.2017.1804

PubMed Abstract | CrossRef Full Text | Google Scholar

Raposo, M., Ribeiro, P., Sério, S., Staiano, A., Ciaramella, A. (Eds.) (2020). Computational intelligence methods for bioinformatics and biostatistics (Cham: Springer International Publishing). Available at: https://link.springer.com/10.1007/978-3-030-34585-3.

Google Scholar

Sasidharan, S. (2012). Modélisation de la croissance des plantes supérieures pour les systèmes de support-vie : modèle métabolique de la feuille de laitue considérant la conversion d’énergie et le métabolisme central du carbone [PhD dissertation]. [Clermont-Ferrand]: Université Blaise Pascal.

Google Scholar

Scattolini, R. (2009). Architectures for distributed and hierarchical model predictive control - a review. J. Process Control 19 (5), 723–731. doi: 10.1016/j.jprocont.2009.02.003

CrossRef Full Text | Google Scholar

Schertl, P., Braun, H. P. (2014). Respiratory electron transfer pathways in plant mitochondria. Front. Plant Sci. 5. doi: 10.3389/fpls.2014.00163

PubMed Abstract | CrossRef Full Text | Google Scholar

Schneider, M., Knuesting, J., Birkholz, O., Heinisch, J. J., Scheibe, R. (2018). Cytosolic GAPDH as a redox-dependent regulator of energy metabolism. BMC Plant Biol. 18 (1), 1–14. doi: 10.1186/s12870-018-1390-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Shameer, S., Ratcliffe, R. G., Sweetlove, L. J. (2019). Leaf energy balance requires mitochondrial respiration and export of chloroplast NADPH in the light. Plant Physiol. 180 (4), 1947–1961. doi: 10.1104/pp.19.00624

PubMed Abstract | CrossRef Full Text | Google Scholar

Soussana, J.-F. X., Graux, A.-I., Tubiello, F. N. (2010). Improving the use of modelling for projections of climate change impacts on crops and pastures. J. Exp. Bot. 61 (8), 2217–2228. doi: 10.1093/jxb/erq100

PubMed Abstract | CrossRef Full Text | Google Scholar

Tan, X. L. J., Cheung, C. Y. M. (2020). A multiphase flux balance model reveals flexibility of central carbon metabolism in guard cells of C3 plants. Plant J. 104 (6), 1648–1656. doi: 10.1111/tpj.15027

PubMed Abstract | CrossRef Full Text | Google Scholar

Taniguchi, M., Miyake, H. (2012). Redox-shuttling between chloroplast and cytosol: Integration of intra-chloroplast and extra-chloroplast metabolism. Curr. Opin. Plant Biol. 15 (3), 252–260. doi: 10.1016/j.pbi.2012.01.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Tatjewski, P. (2017). Offset-free nonlinear model predictive control with state-space process models. Arch. Control Sci. 27 (4), 595–615. doi: 10.1515/acsc-2017-0035

CrossRef Full Text | Google Scholar

Thornley, J. H. M. (2002). Instantaneous canopy photosynthesis: Analytical expressions for sun and shade leaves based on exponential light decay down the canopy and an acclimated non-rectangular hyperbola for leaf photosynthesis. Ann. Bot. 89 (4), 451–458. doi: 10.1093/aob/mcf071

PubMed Abstract | CrossRef Full Text | Google Scholar

Thornley, J., Johnson, I. (1980). Plant and crop modelling: A mathematical approach to plant and crop phisiology (Oxford: The Blackburn Press).

Google Scholar

Tian, X., Wang, P., Huang, D., Chen, S. (2007). Offset-free multistep nonlinear model predictive control under plant-model mismatch. Int. J. Adapt. Control Signal Process. 21, 731–744. doi: 10.1002/acs

CrossRef Full Text | Google Scholar

Vasquez, J., Guerrero, J., Miret, J., Castilla, M., Vicuña, L. (2010). Hierarchical control of intelligent microgrids. IEEE Ind. Electron. Mag. 4 (4), 23–29. doi: 10.1109/MIE.2010.938720

CrossRef Full Text | Google Scholar

Xiao, Y., Chang, T., Song, Q., Wang, S., Tholen, D., Wang, Y., et al. (2017). EPlant for quantitative and predictive plant science research in the big data era–lay the foundation for the future model guided crop breeding, engineering and agronomy. Quant. Biol. 5 (3), 260–271. doi: 10.1007/s40484-017-0110-9

CrossRef Full Text | Google Scholar

Yu, C., Liang, Z., Huang, A. H. C. (1984). Glyoxylate transamination in intact leaf peroxisomes. Plant Physiol. 75 (1), 7–12. doi: 10.1104/pp.75.1.7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, W., Liang, X., Li, K., Dai, P., Li, J., Liang, B., et al. (2021). Metabolomics analysis reveals potential mechanisms of phenolic accumulation in lettuce (Lactuca sativa l.) induced by low nitrogen supply. Plant Physiol. Biochem. 158, 446–453. doi: 10.1016/j.plaphy.2020.11.027

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, J., Wang, J. Z., Hang, T., Li, P. P. (2020). Photosynthetic characteristics and growth performance of lettuce (Lactuca sativa l.) under different Light/Dark cycles in mini plant factories. Photosynthetica 58 (3), 740–747. doi: 10.32615/ps.2020.013

CrossRef Full Text | Google Scholar