TYPE Original Research PUBLISHED 20 October 2022 DOI 10.3389/fpls.2022.970410 Enhancing control systems of OPEN ACCESS higher plant culture chambers EDITED BY Meiyappan Lakshmanan, Bioprocessing Technology Institute via multilevel structural (A*STAR), Singapore REVIEWED BY mechanistic modelling Sanu Shameer, University of Oxford, United Kingdom Lokanand Koduru, Carles Ciurans 1*, Josep M. Guerrero 2, Ivan Martı́nez-Mongue 3, Technology and Research (A*STAR), Singapore Claude G. Dussap 4, Igor Marin de Mas 3 and Francesc Gòdia 1,5 *CORRESPONDENCE 1 Micro-Ecological Life Support System Alternative (MELiSSA) Pilot Plant-Claude Chipaux Laboratory, Carles Ciurans Universitat Autònoma de Barcelona, Barcelona, Spain, 2 Centre for Research on Microgrids (CROM), Carles.Ciurans@uab.cat Aalborg University, Aalborg, Denmark, 3 AAU Energy, Novo Nordisk Foundation Center for Sustainability, Lyngby, Denmark, 4 Institut Pascal, Université Clermont Auvergne, Clermont-Ferrand, SPECIALTY SECTION France, 5 Centre for Space Studies and Research - Universitat Autònoma de Barcelona (CERES-UAB), This article was submitted to Institut d’Estudis Espacials de Catalunya, Universitat Autònoma de Barcelona, Barcelona, Spain Plant Biophysics and Modeling, a section of the journal Frontiers in Plant Science RECEIVED 15 June 2022 Modelling higher plant growth is of strategic interest for modern agriculture as ACCEPTED 20 September 2022 PUBLISHED 20 October 2022 well as for the development of bioregenerative life support systems for space applications, where crop growth is expected to play an essential role. The CITATION Ciurans C, Guerrero JM, Martı´nez- capability of constraint-based metabolic models to cope the diel dynamics of Mongue I, Dussap CG, Marin de Mas I plants growth is integrated into a multilevel modelling approach including mass and Gòdia F (2022) Enhancing control systems of higher plant culture and energy transfer and enzyme kinetics. Lactuca sativa is used as an exemplary chambers via multilevel structural crop to validate, with experimental data, the approach presented as well as to mechanistic modelling. design a novel model-based predictive control strategy embedding metabolic Front. Plant Sci. 13:970410. doi: 10.3389/fpls.2022.970410 information. The proposed modelling strategy predicts with high accuracy the COPYRIGHT dynamics of gas exchange and the distribution of fluxes in the metabolic © 2022 Ciurans, Guerrero, Martı´nez- network whereas the control architecture presented can be useful to manage Mongue, Dussap, Marin de Mas and higher plants chambers and open new ways of merging metabolome and Gòdia. This is an open-access article distributed under the terms of the control algorithms. Creative Commons Attribution License (CC BY). The use, distribution or KEYWORDS reproduction in other forums is permitted, provided the original closed ecological life support systems, higher plant chambers, functional-structural author(s) and the copyright owner(s) modelling, model predictive control, FBA, Melissa are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms. 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 Frontiers in Plant Science 01 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 and agricultural systems on carbon dioxide capture, to mention higher plant growth, surpassing the Farquhar model (Farquhar a few (Soussana et al., 2010; Peng et al., 2020). Besides et al., 1980) to calculate photosynthesis rates based on the responding to the current global climate and demographic enzyme kinetics of principal metabolic pathways. However, challenges, the need for more efficient forms of horticulture to out of the photosynthesis process, many key mechanisms like increase productivities, improve yield, and optimize crop growth respiration, substrate accumulation and management, tissue has also contributed to the generation of mathematical models to morphology, or multi-tissue interactions have not completely support agronomic activities. Still, much progress is required in been mathematically characterised yet. To treat such a complex the field of biological systems modelling and this is especially system, the use of metabolome information has attracted the relevant in the case of higher plants due to the complexity of attention of plant modellers as an alternative to gather multiple their underlying growth mechanisms. Modelling complex biological reactions formalized as a constraint-based metabolic systems like higher plants may be an objective by itself model (Gomes et al., 2015). It should be highlighted that as before-mentioned, but they may also lead to other constraint-based metabolic models have been demonstrated to interesting applications like the development of model-based successfully address the plant diel cycle with a light phase with control methodologies. Particularly, space research and the resource accumulation and a dark phase with resource depletion development of bio-regenerative life support systems (BLSSs), (Cheung et al., 2014), a critical phenomenon in higher plant which are the set of technologies designed to guarantee life in metabolism very difficult to deal through first-principle long-term crewed missions (Eckart, 1995), have exploited the approaches thus far. Several efforts have been recently placed use of controllers based on first principles in opposition to on integrating available mechanistic information and omics data surrogated and reduced order models, to improve both the in a common multiscale modelling framework that could management of missions under operation as well as to potentially be used by the plant computational biology improve the design of future missions (Fulget et al., 1999). research community to feed data in a single converging One of the main actors on BLSS research is the Micro- platform (Marshall-Colon et al., 2017; Xiao et al., 2017). Ecological Life Support System Alternative (MELiSSA), a In this study, the modelling of higher plants is approached European Space agency (ESA) BLSS program, which is through the design of a multilevel organization of the devoted to developing technologies accompanied by a mechanistic processes that take place during crop growth. To modelling framework to support their research and present the results, L. sativa has been used as an exemplary development activities. One remaining task, which is in turn higher plant. The higher level in the hierarchy copes with the one of the most relevant consensuses of the development of mechanistic phenomena corresponding to a higher BLSSs shared by the major space agencies, is the importance of characteristic length (i.e., crop chamber scale), whereas the developing technologies and mathematical models to grow lower the level the smaller the characteristic length of the plants on space, which are expected to be the major source of modelled phenomena (i.e., enzyme rate). Information follows a edible biomass in BLSSs (Gitelson and Lisovsky, 2002; Poughon top-bottom flow, and it is eventually used to calculate the et al., 2009; Dong et al., 2017). metabolic flux distribution by applying a flux balance analysis Even though current models can cope with the evolution of (FBA). This multilevel modelling approach is firstly validated biomass, the compounds involved in photosynthesis and with experimental data and secondly integrated in a model- respiration (O2 and CO2), and the nutrient uptake by the based predictive control, representing, to the best of the authors, roots, most of the phenomena involved in plant growth are the first attempt to incorporate cell metabolism in an advanced not addressed given their complexity and lack of knowledge. control strategy. This complexity is associated with several factors; higher plants Overall, the modelling and control methodology presented are multicellular, compartmentalized organisms undergoing in this study may pave the way for a more efficient and strong metabolic changes associated with the cyclic switch sustainable agriculture either for intensive cultivation systems between light and dark phases of the day. They are also or as a part of BLSSs in space exploration. characterised by having a complex substrate partitioning strategy with different parts being coordinated to uptake and distribute specific compounds. This has contributed to the 2 Model proposal preference of empirical at the expense of mechanistic models due to the usually satisfactory information provided by the 2.1 Multilevel mechanistic model former, especially in nominal conditions (Boscheri et al., 2012; Amitrano et al., 2020). Empirical models cover a limited range of The model developed here is organized following a operating conditions though; thus, the scope of their use is multilevel approach, considering the different levels of the narrow and cannot contemplate all scenarios that plant culture plant, from canopy to metabolic level, and uses the output of may undergo. As an alternative, mechanistic approaches have the higher levels as the input to the lower levels. In this section, also been deployed to understand the first principles behind the models used in the different levels are explained. A graphical Frontiers in Plant Science 02 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 description of the model organization is presented in Figure 1 extinction coefficient k: and detailed in the following paragraphs. For the sake of clarity, model parameters are described in Table 1 and constant Is = k · fs · Iu + k · (1 − fs ) · Iu · e(−k·LAI) (2) parameters in Table 2. Id = k · (1 − fs ) · Iu e(−k·LAI) (3) 2.1.1 Level 3: Modelling canopy growth Notice in (2) that the parts of the canopy under direct When dealing with whole-leaf or canopy modelling, it is irradiance also receive diffuse irradiance. The LAI term should necessary to consider the effect of shading among leaves, which be differentiated between the fractions exposed to direct and is not accounted in single-leaf models (Poulet, 2018). The diffuse light sources (LAIs and LAId respectively), as suggested by common practice is to use the leaf area index (LAI) as an Thornley (2002): indicative parameter of the leaf area density over ground LAIs = (1 − e−kLAI )=k (4) surface. The photon flux density inside the canopy I declines along the canopy exponentially and is a function of the leaf area index: LAId = LAI − LAIs (5) The derivative of LAIs and LAId should be obtained and used I = Iu e(−kLAI) (1) to integrate Is and Id to calculate the overall irradiance received Parameter Iu represents the photon flux density at the top of by the canopy: the canopy, and k represents the extinction coefficient. Extending (1), leaves can receive direct photon flux density (Is) dLAIs = e−kLAI dLAI (6) or diffuse photon flux density (Id) as stated by Thornley (2002), which are expressed in terms of mmol m-2 leaf s-1 by using the dLAId = ( 1 − e−kLAI ) dLAI (7) FIGURE 1 Graphical scheme of the multilevel mechanistic model approach used for L. sativa prediction growth. Frontiers in Plant Science 03 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 TABLE 1 Model parameters. Symbol Description Units Morphological module LA Leaf area m2 leaves L Leaf length m leaves Mx Dry biomass g Irradiance module I Irradiance moles m-2 ground s-1 Iu Irradiance at LAI = 0 moles m-2 ground s-1 Is,g Direct irradiance over ground surface moles m-2 ground s-1 Id,g Diffuse irradiance over ground surface moles m-2 ground s-1 Is Direct irradiance over leaf surface moles m-2 leaves s-1 Id Diffuse irradiance over leaf surface moles m-2 leaves s-1 LAI Leaf area index m2 leaves m-2 ground LAIs LAI exposed to direct irradiation m2 leaves m-2 ground LAId LAI exposed to diffuse irradiation m2 leaves m-2 ground Energy balance Tl Leaf temperature K kt Heat transfer coefficient m s-1 Ehs Direct irradiance energy J s-1 Ehd Diffuse irradiance energy J s-1 Er Radiation energy J s-1 Econv Convection energy J s-1 Etr Transpiration energy J s-1 Gas exchange ExCO2 CO2 exchange rate moles m-2 leaves s-1 ExO2 O2 exchange rate moles m-2 leaves s-1 ExH2O H2O exchange (transpiration) rate moles m-2 leaves s-1 Gz Conductance compound z moles m-2 leaves s-1 Pbz Bulk partial pressure compound z Pa Plz Leaf partial pressure compound z Pa Biochemical module FLETC Light electron transport chain rate moles m-2 leaves s-1 max FLETC Maximum light electron transport chain rate moles m-2 leaves s-1 J Ribulose 1,5-biphosphate regeneration moles m-2 leaves s-1 Г Carbon dioxide compensation point moles m-3 Vc Carboxylation rate moles m-2 leaves s-1 Vcmax Maximum carboxylation rate moles m-2 leaves s-1 Vo Oxygenation rate moles m-2 leaves s-1 Ci Carbon dioxide leaf concentration moles m-3 Oi Oxygen leaf concentration moles m-3 Pg Gross photosynthesis rate moles m-2 leaves s-1 Pn Net photosynthesis rate moles m-2 leaves s-1 Boundary layer z gBL Boundary layer conductance of compound z moles m-2 leavess-1 gsz Stomatal conductance of compound z moles m-2 leaves s-1 d Boundary layer thickness m Tbl Average leaf-bulk temperature K Dz Diffusion coefficient of compound z m2 s-1 vbulk Bulk velocity m s-1 vfree Free velocity m s-1 rl Leaf air density kg m-3 Frontiers in Plant Science 04 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 TABLE 2 Model constants. Symbol Description Value Units Reference Irradiance module K Extinction coefficient 0.5 MELiSSA Pilot Plant fs Direct irradiance fraction 0.7 Thornley and Johnson (1980) BCmol C-mole molecular weight 27 g mol-C-1 MELiSSA Pilot Plant Energy balance Na Avogadro number 6.02·1023 Pa c Light velocity 3·109 Pa h Planck constant 6.63·10-34 m2 kg s-1 ϵ Leaf emissivity 0.97 Poulet et al. (2020) s Stefan–Boltzmann constant 5.67·10-8 J s-1 K-4 Poulet et al. (2020) R Ideal gas constant 8.314 m3 Pa K-1 mol-1 Cp Molar air-specific heat capacity 29.3 J mol-1 K-1 Poulet et al. (2020) 4 lmol Water latent heat of vaporization 4.0788·10 J mol-1 Poulet et al. (2020) Biochemical module q Convexity coefficient 0.8 Farquhar et al. (1980) f Energy loss for LETC 0.045 Nikolov et al. (1995) Kc Carboxylation half-saturation constant 460 mbar Farquhar et al. (1980) Ko Oxygenation half-saturation constant 330 mbar Farquhar et al. (1980) max25 FLETC max FLETC at 25°C 100 mmol m-2 s-1 Nikolov et al. (1995) Vcmax25 Vcmax at 25°C 31.31 mmol m-2 s-1 Nikolov et al. (1995) -1 -1 E Activation energy of reaction 81,993 J mol K Nikolov et al. (1995) S Entropy 711.36 J mol-1 K-1 Nikolov et al. (1995) H’ Energy of deactivation 219,814 J mol-1 K-1 Nikolov et al. (1995) Mc C-molar molecular weight 27 g mol-C-1 MELiSSA Pilot Plant DM Dry biomass fraction 0.045 g/g MELiSSA Pilot Plant Boundary layer h Air kinematic viscosity 1.8·10-5 m2 s-1 Poulet et al. (2020) a Leaf angle in relation to the vertical axis 0.1 ° Poulet et al. (2020) g Gravity force 9.8 m s-2 rb Bulk air density 1.186 kg m-3 Z LAI 1 Is,l = Is dLAIs (8) Ehd = Id · Na · c · h · LAId (11) 0 l600nm Z Na represents the Avogadro number, c represents the LAI Id,l = Id dLAId (9) velocity of light, and g is the fraction of wavelength l that 0 compose the light directly irradiating the canopy. The total The way light irradiates exposed and shadowed leaves energy irradiated to the leaves is the summation of both strongly affects the energy balance in the leaf surface, with the Equations (10) and (11). The radiation energy emitted by the shadowed leaves irradiated by diffused light mainly of a plants (Er), the energy lost by convection (Econv), and the energy wavelength of 600 nm corresponding to the green-colour lost by transpiration (Etr) are determined by the following spectrum of transmitted light. Energy received by irradiance equations: contains the direct (Ehs) and diffuse (Ehd) terms, which are   Er = ϵs Tleaf 4 − Tb4 (12) expressed as follows: i=lmax gi Pb (T − Tb ) Ehs = Is · Na · c · h · LAIs o i=l li (10) Econv = Cp kt RTb leaf (13) min Frontiers in Plant Science 05 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 Etr = lmol ExH2 O (14) power. Thus, it is necessary to define an expression for FLETC, which can be approximated by a quadratic equation (Nikolov In (12), ϵ and s represent the leaf emissivity and the Stefan– et al., 1995): Boltzmann constant, respectively. In (13), Cp and kt represent the qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi molar specific heat capacity at constant pressure and 298.15 K and max (FLETC + jIu ) − max (FLETC + jIu )2 − 4jqIu FLETC max the heat transfer coefficient, respectively, the latter being a function FLETC = (19) 2q of the diffusion coefficient and the boundary layer thickness, as follows (see Supplementary S1 for details on Dt calculation and for (1 − f ) an extended description of the boundary layer model): j= (20) 2 Dt In (19), the maximum FLETC is represented by FLETC max (see kt = (15) d Supplementary S1 for calculation details). The photosynthetic In (14), lmol is the water latent heat of vaporization and photon flux density used is represented by Iu and q is a convexity ExH2O the transpiration rate defined in the following section. coefficient. The efficiency of energy conversion is represented by Finally, Tleaf is the leaf surface temperature, Tb is the bulk j, which is a function of the fraction of absorbed photon flux temperature, and R is the ideal gas constant. unavailable for photosynthesis (f). As previously demonstrated, it is necessary to consider the direct and diffuse irradiances when 2.1.2 Level 2: Modelling gas exchange rates considering the whole canopy. Therefore, combining with (4) This level is dedicated to calculating the uptake and release and (5): rates between the atmosphere and the leaves concerning the Z LAI FLETC (Is ) dLAIs exchange gases (ExCO2 , ExO2 , and ExH2O ). The approach to FLETC(Is ,Id ) = 0 2 modelling gas exchanges between the leaves and the atmosphere Z LAI follows Fick’s law, with the concentration gradient being the FLETC (Id ) dLAId + (21) driver of the molecular transport: 0 2 GCO2 CO2 Considering that two electrons are necessary per molecule of ExCO2 = (Pb − PlCO2 ) (16) NADPH generated, the light electron transport chain rate Pb resulting from (21) is divided by 2. This is necessary to provide information to the FBA with consistent units GO2 O2 ExO2 = (P − PlO2 ) (17) considering the stoichiometry matrix used (see S2 for a full list Pb b of reactions and stoichiometry matrix). As mentioned before, in Farquhar et al.’s (1980) approach, it is GH2 O H2 O ExH2 O = (Pb − PlH2 O ) (18) necessary to convert FLETC defined in (21) to a flow of ribulose 1,5- Pb biphosphate (Ru5P) regeneration (J), through the following In (16)–(18), the atmospheric partial pressure is calculated expression derived from the stoichiometry of the light electron assuming gases behave following the general gas equation, transport chain and the Calvin cycle (Farquhar et al., 1980): whereas the conductance for the different gases (GCO2 , GO2 , FLETC and GH2O ) and the internal (i.e., leaf) partial pressure (PlCO2 , PlO2 , J= (22) 2(2 + 2 ∅ ) PlH2 O ) as well as the leaf area (LA) are calculated according to Poulet et al. (2018) and explained in Supplementary S1. The RuBisCo carboxylation (Vc) kinetics is of Michaelis– Menten type and is a function of the leaf internal oxygen (Oi) 2.1.3 Level 1: Modelling enzyme kinetics and carbon dioxide (Ci) concentrations: The gross photosynthesis rate (Pg) is calculated using the Farquhar model (Farquhar et al., 1980), which has been widely (Ci − G ) Vc = Vcmax (23) used to model photosynthesis phenomena (Harley et al., 1992; Ci + Kc (1 + KOoi ) Morgan and Robles, 2002; Arnold and Nikoloski, 2011). In this model, Pg is determined through finding the limiting rate of the In (23), Vcmax represents the maximum carboxylation photosynthesis, which is caused by either the regeneration of velocity of RuBisCo, K c and K o are the Michaelis– ribulose 1,5-biphosphate (J), substrate of the ribulose-1,5- Menten half-saturation constants for the carboxylation and biphosphate carboxylase (RuBisCo), or the RuBisCO oxygenation activities of RuBisCo, respectively, and Г carboxylation rate itself (Vc). represents the carbon dioxide compensation point. Details On the one hand, the regeneration of ribulose 1,5- about the calculation of the internal carbon and oxygen biphosphate depends on the potential rate of the light electron concentrations are found in Supplementary S1, based on the transport chain (FLETC) and its capacity to generate reducing boundary layer approach defined by Poulet et al. (2018; 2020). Frontiers in Plant Science 06 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 Vc can be directly used to feed the FBA as an upper bound 2.1.4 Level 0: Stoichiometry matrix (see Figure 1). In this level, a simplified network model of the Finally, the gross photosynthesis rate (Pg) is determined by photosynthetic leaves’ metabolism of L. sativa is described. finding the minimum between the Ru5P regeneration rate (J), The stoichiometry matrix is based on the work of Sasidharan the RuBisCo carboxylation rate (Vc), and the gas exchange rate (2012), which contains the distinguishing characteristics of L. (ExCO2 ) and the net photosynthesis rate (Pn) is determined by sativa, such as the reduced starch content to store carbon and the retrieving the RuBisCo oxygenation (Vo) to Pg as follows: definition of the elemental composition that makes up the macromolecules of the biomass. This model though has been Pg = min(Vc , ExCO2 , J) · LA (24) extended to include relevant reactions like the pentose phosphate pathway or the photorespiration cycle, originally Ol Vomax Ko missing. In the model used, the cellular organelles are Vo = Vc · · Cl (25) Vcmax Kc 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 Pn = Pg − Vo · LA (26) dark phases of the day account for separate pools of metabolites Notice that Equation (24) can be formalized due to the and organelle compartments. Here, only those metabolites that conversion applied in Equation (22) considering that J have been reported to be accumulated in one phase and represents the regeneration rate of Ru5P and Vc represents its consumed in the other are connected by exchange reactions carboxylation rate. Parameter A indicates the surface of the crop that simulate the transference of nutrients between phases. For growing area. In this study, the use of a metabolic matrix makes example, sugars that are synthesised in the light phase can be the conversion from light electron transport flux (FLETC) to used in the dark phase due to the addition of exchange reactions RuBP regeneration (J) unnecessary, because this information is among day phases. In Figure 2, the metabolic model is presented already included in the stoichiometric matrix. Similarly, the in a simplified way including the cellular compartmentalization, discontinuity introduced by Farquhar et al. (1980) in (24) can the pathways involved, and the connections between them as be prevented with an FBA formulation as addressed in the well as the exchange reactions with the atmosphere and between following section. day periods. 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 (CO2, O2, H2O, light photons) and through the roots (HNO3 and H2SO4). PG represents 2-phospho-glycolate. Frontiers in Plant Science 07 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 The diel model is composed by different organelles including represents the array of fluxes for each of the reactions of the chloroplasts, mitochondria, peroxisomes, and cytoplasm. The model. The lower and upper bounds (lb and ub) are fixed only proposed model also features plants’ simultaneity of metabolic for those fluxes indicated in Table 3. Matrix C contains pathways (i.e., glycolysis in chloroplast and cytoplasm or the information regarding the flux ratios specified in Table 3, and folate metabolism in chloroplasts, cytoplasm, and mitochondria) matrix F represents the resulting flux. The inequality constraints as well as the aforementioned coordination between day represented in Table 3 are generated in levels 3 and 2 of the phases. A summary of the metabolic model is found in model and are used to feed the FBA. ATPmaintenance and NAD(P) Supplementary S3, and the model in SBML format is accessed Hmaintenance includes the reactions that contribute to the in Supplementary S4. consumption of ATP and the reducing agent for respiration The metabolic network model is mathematically formalized purposes. As suggested in Cheung et al. (2013), this can be as a constraint-based model, and the fluxes are calculated by achieved including generic ATPase and NAD(P)H oxidase applying a flux balance analysis (Cheung et al., 2014): reactions. Finally, the enzyme rates included are those related to constrain the reducing power supply in the cytoplasm. At max n i g h t , p l a s t i d ic N A D P - m a l a t e d e h y d r o g e n a s e a n d (vbiomass ½d) v glyceraldehyde 3-phosphate dehydrogenase are downregulated (27) lb ≤ Sv ≤ ub (Miginiac-Maslow and Lancelin, 2002; Mekhalfi et al., 2014). subject : Finally, minimums and maximums for a set of reactions are Cv = F defined given the presence of thermodynamically infeasible In the FBA formulation in (27), the objective function refers loops when no restrictions are applied. These reactions include only to the daily biomass production. It is well known for other the following: PYK, pyruvate kinase; PGM: phosphoglycerate crop species like Arabidopsis thaliana that biomass growth also mutase; ENO: enolase; EB1: inorganic pyrophosphatase; EB2: takes place during the dark phases of the day (Gomes et al., inorganic pyrophosphatase; ACS: acetyl-CoA synthetase; 2015). However, night metabolism of non-starchy crops like Ser_bio_cl: phosphorylated serine pathway; GOGAT: lettuce is still not clear so in the current approach maintenance- glutamate synthase; Prot32: 3-mercaptopyruvate associated reactions are limited to night metabolism whereas sulfurtransferase/cytoplasmic aspartate aminotransferase; metabolism associated with light periods of the day concentrates, OASTL: cysteine synthase; GS: glutamine synthetase. The FBA on top of maintenance, also biosynthetic reactions. Letter υ presented is implemented both in Matlab ® 2021 using the TABLE 3 Flux balance analysis equality, inequality, and flow ratio restrictions. [d] and [n] indicate day and night phase period respectively. Inequality constraints Type Description nExCO2 lb Gas exchange nExO2 ub Gas exchange nExH2O lb Gas exchange nFLETC ub Light ETC nVc ub Carboxylation nVo lb Oxygenation Flux ratios Value nExO2[d] : nExCO2[d] 1.22:-1 Photosynthesis rate (MELiSSA Pilot Plant) nATPmaintenance[d/n]: nNADPHmaintenance[d/n] 3:1 Maintenance (Cheung et al., 2014) nExCO2[d] : nExCO2[n] -1:0.25 Respiration (MELiSSA Pilot Plant) nATPmaintenance[d]:nATPmaintenance[n] 1 Respiration (MELiSSA Pilot Plant); (Liu and van Iersel, 2021) nNADPHmaintenance[d]:nNADPHmaintenance[d] 1 Respiration (MELiSSA Pilot Plant); (Liu and van Iersel, 2021) nOPPP[c,m,cl][d/n]+nICDH[c,m][d/n]+nME[c,m][d/n]:nNADPHmaintenance[d/n] 1:1 (Corpas and Barroso, 2014) nSGT[p][d/n] : nGT[p][d/n] 1:1 (Yu et al., 1984) Enzyme rates Value nGAPN[c] ub: 0.33 (Shameer et al., 2019) nGAPDH[c] lb: -93 (Shameer et al., 2019) nMDH[c] lb: -0.75 (Shameer et al., 2019) *NADPH oxidation and ATP hydrolysis associated with maintenance reactions include plastidic, cytoplasmic, and mitochondrial locations. Subscripts c, m, cl, and p indicate cytoplasmic, mitochondrial, plastidic, and peroxisomal location. lb and ub represent lower and upper bounds, respectively. List of enzyme abbreviations: GAPN, cytosolic non-phosphorylating NADP-glyceraldehyde-3-phosphate dehydrogenase; GAPDH, cytosolic glyceraldehyde 3-phosphate dehydrogenase; MDH, cytosolic malate dehydrogenase; OPPP, glucose-6-phosphate dehydrogenase and 6-phosphogluconate dehydrogenase; ICDH, isocitrate dehydrogenase; MDH, malate dehydrogenase; SGT, serine- glyoxylate aminotransferase; GT, glutamate-glyoxylate aminotransferase. Frontiers in Plant Science 08 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 Cobra Toolbox (Becker et al., 2007) and in python 3.0 using estimation of the controlled variables to the setpoint at steady Cobrapy (Ebrahim et al., 2013) and can be found in state. This output becomes the input of the following control step Supplementary S5. in the hierarchy, which is a model-based predictive control (MPC). MPC uses a discretized model of the process to be 2.1.5 Dynamic model controlled and aims at finding the sequence of control The evolution of the different states of interest through time commands along a prediction horizon that brings the predicted for either dry biomass, leaf temperature, or gas compositions can controlled variables the closest to the reference given a set of hard be obtained integrating their rates of generation or consumption and soft constraints. MPC works based on a rolling-horizon over time. For dry biomass, this can be done straightforward approach, which essentially solves the minimization problem from Pn, considering that all carbon molecules captured by the along the defined prediction horizon but sends solely the plant are fixed into structural biomass: control command corresponding to the first step. This process Z is repeated every time the controller is executed (Pannocchia and t= tf Mx = Pn · BCmol · dt (28) Bemporad, 2007). In Figure 3, a schematic representation of the t=t0 advanced control architecture with details on the communication For the evolution of leaf temperature, it is necessary to solve among its levels is presented. In the presented study, the an energy balance between the leaf temperature and the controlled variable is the oxygen concentration in the chamber environment and to convert energy units to temperature whereas the manipulated variable is the gas flow. degrees: The SSTO solves the following system of equations using the output of FBA, considering oxygen exchange rate (vExO2 ) as the Z t= tf (Ehs + Ehd )A − (Er + Econv + Etr )LA generation rate in a steady-state mass balance: Tl = Mx · dt (29) t=t0 Cp DM 2 − xs ) + vExO2 = 0 us (Oin (31) Finally, the oxygen concentration can be obtained by solving a mass balance within the growing crop chamber: xs + d = ys (32) Z t= tf u· (Oin − O2 ) + vExO2 2 O2 = · dt (30) ys = yref (33) t=t0 V In (30), gas flow is represented by u, oxygen concentration in s.t. the input flow by Oin2 , and the chamber volume is V. yL ≤ ys ≤ yU (34) uL ≤ us ≤ uU (35) 2.2 Integration of the multilevel mechanistic model with advanced As indicated in Figure 3, the output of the SSTO is provided control architectures as a reference to the MPC, which in this case is the concentration of oxygen in the chamber (xs) and the external gas flow (us) both Once having defined the modelling strategy for higher crops, at steady state. The internal model at steady state is defined in the second scope of this study is to integrate the use of metabolic (3), with the generation term (vExO2 ) being the output of the models into an advanced control strategy. In different complex FBA, Oin2 the input oxygen concentration, and the internal model systems ranging from microgrids (Vasquez et al., 2010), life prediction defined by xs. A disturbance (d) is incorporated in support (Ciurans et al., 2021), or water distribution systems (32) to take into account any possible plant-model mismatch or (Ocampo-Martinez et al., 2012) to chemical plants (Scattolini, a measured perturbation and must be taken into account for the 2009; Marchetti et al., 2014), the use of advanced control new reference generation. Hence, the final prediction value (ys) architectures has been proven to be efficient in terms of optimal at steady state considering the presence of any given disturbance management and control. Advanced control architectures are will match the process measurement guaranteeing offset-free characterized by hierarchically distributing management and control (Pannocchia and Bomperad, 2007). The expression to control functions in different levels. In this study, an adaptation obtain the disturbance is defined in the MPC development of a common control architecture used in process plants is hereafter. The technical upper and lower bounds of controlled adapted to control oxygen in a crop-growing chamber. The top and manipulated variables are summarised in (34)–(35) and layer called steady-state target optimization (SSTO) aims at thus do not need to be defined in the MPC. It is then possible to finding reference values for the controlled and manipulated violate constraints on controlled and measured variables during variables given a specific setpoint through solving a mass transition states but not at steady state (Marchetti et al., 2014). balance problem at steady state. The output of SSTO are the The MPC solves a rolling-horizon non-linear optimization controlled and manipulated variables that give the closest problem, taking the output of the SSTO as the reference to track: Frontiers in Plant Science 09 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 FIGURE 3 Scheme of the advanced control architecture proposed to integrate metabolic models. i=Np jx^ ðijkÞ − xs j2 i=Nc ju(ijk) − u j2 minimum values for both controlled and manipulated variables J = m in l1 o + l2 o s (defined by Dx and Du) and are subject to a scaling factor (l). U i=k Dx i=k Du Prediction and control horizons are represented by Np and Nc, i=Nc ju(ijk) − u(i − 1jk)j2 respectively. Constraints in (37) reflect the dynamics of the growing + l3 o (36) i=k Du crop chamber with Ts as the sample time of the MPC. The internal model is initialized with the current process measurement (yp). In s.t. (39), the disturbance is integrated to the internal model prediction    T T as similarly done in the SSTO in (32). This disturbance is estimated x^ (i + 1jk) = x^ (ijk) + U(ik) Oin ^ (ijk) s + vExO2 s 2 −x (37) V V at each sampling time and is defined as the difference between the process measurement (yp) and the first step prediction of the y^ (ijk) = x^ (ijk) + dðijk) (38) previous MPC execution (^ y (2jk − 1)) 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 dðik) = yp − x^ (2k − 1) (39) 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 dði + 1jk) = dðijk) (40) the next SSTO and MPC execution. This control strategy is 2 3 implemented in Matlab® 2021 using the Optimization Toolbox U(i ∣ k) for solving non-linear programming problems and Cobra Toolbox 6 7 U=6 4 ⋮ 7 5 (41) 2022 for the FBA resolution (Becker et al., 2007). U(Nc ∣ k) The cost function in (36) includes penalization terms to the 2.3 Simulation scenarios deviation of the internal prediction (^ x ) and the control command (u) from the concentration (xs) and the gas flow (us) references Two simulation packages are presented: first, the results of generated in the SSTO. A third penalization term to the rate of the multilevel model presented in Section 2.1 and their change of the manipulated variable is also included in (36) aimed to validation with experimental data, and second, a dynamic adjust the speed of the controller response. All three penalization simulation presenting the response of the control architecture terms are normalized using the range of possible maximum and presented in Section 2.2 under different perturbations. Frontiers in Plant Science 10 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 2.3.1 Simulation conditions for multilevel 3 Results model validation CO2-response curves were generated using the fixed light This section is divided into two parts: the first dedicated to intensity indicated in Table 4 and the following range of internal present the output of the multilevel model and its validation with CO2 values in mmol CO2/mol: 100, 225, 300, 450, 600, 850, 1,000, experimental data and the second related to the integration of 1,100. Light-response curves were generated using the fixed that modelling approach containing metabolic information into internal carbon dioxide concentration indicated in Table 4 and an advanced control architecture. the following range of light intensity values in mmol/m2/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 3.1 Validating photosynthesis rates to retrieve the inequality constraints described in Table 3 and thus the FBA can be resolved. Results are represented in Section The results of the model introduced in Section 2.1 are 3.1. The distribution of fluxes is analysed for both light and dark represented in Figure 4, where modelled and experimental metabolism using atmospheric conditions for CO2 which are results for mature leaves (28 days after transplanting) 400 ppm and a light intensity of 400 μmole/m2/s. Results are are compared. graphically represented in Section 3.2. Figure 4A represents the RuBisCo saturation curve showing a fast rate of change in the smaller range of internal carbon 2.3.2 Advanced control architecture dioxide concentrations, a pattern which is reproduced in configuration Figure 4B which shows the light electron transport chain In Table 5, the controller specifications and parameter values saturation curve. The maximum net photosynthesis rate are indicated. On top of the control objective defined in Section achieved in the light-response curve is lower than the 2.2, atmospheric CO2 is controlled at 800 ppm with external maximum achieved in the CO2-response curve because in the addition of pure CO2. former, the internal CO2 concentration used for the simulation Worthy of note is that the scaling factors are used to promote is 400 mmol CO2/mol reaching the expected photosynthesis rates the control of the system close to the setpoint but at the expense of if compared to Figure 4A. Finally, the output of the model is having a more aggressive control. Prediction and control horizons comparable to the reported experimental results by Zhou et al. are important tuning parameters of the controller increasing its (2020) under the same operating conditions. Overall, the error sensitivity but also the computational cost of the calculation. Finally, observed in the modelled results in relation to the experimental the sampling time also affects the control performance. For slow values is higher at low internal carbon concentrations and at low systems like a plant cultivar, sample times should not be too short light intensities. 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 3.2 Distribution of metabolic fluxes using trade-off between the expected performance and the controller and a metabolic diel model system capabilities. The metabolic-based control architecture represented in Figure 3 is tested in a 24-h dynamic simulation In this section, the flux distribution of day and night using the operating conditions listed in Table 4. metabolites obtained after the resolution of the FBA TABLE 4 List of operating parameters used in FBA simulation. Parameter Value Units Leaf area (LA) 25 m2 leaf Leaf area index (LAI) 5 m leaf/m2 ground 2 Growing area (A) 5 m2 ground Chamber height (H) 1 m Bulk temperature (Tb) 25 °C Bulk pressure (Pb) 101,300 Pa Relative humidity (RH) 70 % Light intensity (Iu) 800 mmol/m2/s Internal CO2 concentration (Ci) 400 mmol/mol Forced velocity (vforced) 0.3 m s-1 Frontiers in Plant Science 11 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 TABLE 5 Controller specifications including SSTO and MPC algorithms. Parameter Value Units Input oxygen concentration (Oin 2 ) 0 % Oxygen setpoint (yref) 21 % Carbon dioxide setpoint 800 ppm Lower and upper bounds oxygen concentration (yL – yU) 18-24 % Lower and upper bound flow ((uL – uU)) 0-1 m3/h Scaling factor (l1,l2,l3) 10,1,1 Prediction horizon (Np) 4 Control horizon (Nc) 3 Sample time (Ts) s introduced in (27) are presented. The parameters used for the As previously stated, higher plants store carbon during the simulation are those provided in Table 3. light phase of the day, to be used during night respiration and fuel maintenance processes. For non-starchy crops, even 3.2.1 Day flux distribution though they can still generate starch, most of the carbon The day flux distribution is presented in Figure 5. The fixed during the light phase is stored as soluble sugars or central carbon metabolism of plants in light conditions is well organic acids. In this study, the sugar molecules stored and represented in this model with the main fluxes located in the mobilized between light and dark periods of the day have not Calvin cycle-associated reactions. The results indicate that, as been restricted and, as Figure 5 suggests, sucrose, fructose, has been extensively studied and published (Michelet et al., 2013; citrate, malate, and fumarate are the metabolites used for Tan and Cheung, 2020), the flow through the Calvin cycle carbon exchange. This modelled result fits well with the generates triose phosphate (g3p) from ribulose 1,5- reported experimental concentrations of sugars in lettuce at biphosphate (RuBP), consuming part of the reducing power harvest, with glucose, sucrose, and fructose being the main and ATP molecules synthesized during the light electron carbohydrates found for carbon exchange between phases of transport chain. Triose phosphate is used to feed the rest of the day (Chen et al., 2019) and also predicting with accuracy the Calvin cycle machinery aimed at restoring the ribulose 1,5- the role of malate accumulation during the light phase of the biphosphate while it is also partially used to generate day in vacuoles for its use in the dark (Lee et al., 2021). Not all photosynthetic end products (McClain and Sharkey, 2019). carbon compounds mobilized from light to night metabolism A B 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 mmol/m2/s is set; in (B), a fixed internal carbon concentration of 400 mmol/m2/s is set. Experimental data are based on Zhou et al. (2020). Frontiers in Plant Science 12 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 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: FLETC, 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; Vc, RuBisCo carboxylation; Vo, 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. are consumed during the latter. Therefore, some carbon One of the critical phenomena of plant photosynthetic cell intermediates need to be exported from night to light metabolism is the coordination of photosynthesis and metabolism too. Specifically, citrate imported from dark respiration, which essentially determines how and where periods is used in the light phase of the day to generate energy carrier molecules (ATP and NAD(P)H) are produced. oxoglutarate which is important for the nitrogen assimilation Most of the ATP and NAD(P)H used for catabolic reactions are mechanism and for the synthesis of nitrogen-rich amino acids. produced in the light electron transport chain in chloroplasts for The citrate cycle and its interactions with amino acid amino acid and lipid production. Part of the ATP synthesized in biosynthesis are well covered by the presented model, both chloroplasts is exported to the cytosol through the 3PG-G3P suffering a flux reduction when nitrate uptake is limited shuttle, satisfying the ATP demand together with ATP exported (Morcuende et al., 1998; Weiwei Zhou et al., 2021). from mitochondria (Gakière et al., 2018; Shameer et al., 2019). Frontiers in Plant Science 13 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 The resulting metabolic network shows the mechanisms of light photoperiods. This is mainly because most TCA redox power balancing in the different organelles of intermediates are dedicated to anabolic reactions and pyruvate photosynthetic leaves enabled by the metabolite shuttles dehydrogenase (PDH) is photo-inhibited (Schertl and represented by the malate/oxalacetate and the triose Braun, 2014). phosphate/3-phosphoglycerate shuttle in Figure 5 (and the glutamate/2-oxoglutarate and the malate/aspartate shuttle, not 3.2.2 Night flux distribution represented) (Taniguchi and Miyake, 2012). Around half the During the dark phase of the day in Figure 6, ATP can only NAD(P)H generated in mitochondria comes from the glycine be produced in mitochondria with NAD(P)H being the electron decarboxylation, which in turn generates the serine used in the donor and thus completely modifying the flow distribution serine-glyoxylate aminotransferase (SGT) in the peroxisome. It within the cell. Autotrophic organisms like plants use the has been reported that malate dehydrogenase seems to regulate sugars generated and stored during the day to feed TCA, the reducing power in mitochondria based on the reduction state which is cyclic in dark conditions. In lettuce and other non- of the cells, removing and restoring NAD(P)H at low- and high- starchy vegetables, instead of mobilizing starch, glucose, sucrose, light conditions responding to changes in the photosynthesis or fructose is broken down to pyruvate at very similar rate (Bykova et al., 2014; Schertl and Braun, 2014). This is proportions (Chen et al., 2019) even though in this study the validated in the presented fluxome, where mitochondrial malate proportion of soluble sugar utilization has not been constrained. dehydrogenase (MDHm) removes excess NAD(P)H produced The pyruvate produced in the glycolysis is then transferred to the through the glycine cleavage system (GCS, also known as glycine mitochondria to regenerate the reducing power needed to fuel decarboxylase system), the main contributor of redox power in the mitochondrial respiration. In dark conditions, the exchange the mitochondria. The metabolic flux distribution represented in rates are completely opposite to those observed in light Figure 5 also shows that the TCA cycle is not complete during conditions, with carbon dioxide and water being released and 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, a-oxoglutarate dehydrogenase; SUCOASm, succinyl-CoA synthetase; SUCD1m, succinate dehydrogenase complex; MDHm, malate dehydrogenase; MP1/2, mitochondrial phosphorylation 1/2. Frontiers in Plant Science 14 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 oxygen consumed. In a way, respiration and photosynthesis are metabolism and is represented in Figure 7 considering an opposite processes, but complementary as demonstrated during atmospheric concentration of CO2 of 1,000 ppm in order not crop growth. The governing rule of night metabolism is to limit by carbon substrate and a nominal light intensity of 400 represented in Table 3 in the ratio of consumed CO2 during μmol/m2/s. the day phase over produced CO2 due to night respiration Considering the light metabolism sensitivity analysis (ExCO2[d]:ExCO2[n]) which sets the night carbon conversion. (Figure 7A), the light phase of day-time mitochondrial respiration is never active. In the predicted fluxome for all 3.2.3 Sensitivity analysis ranges of flux ratios applied, the ATP production comes FBA is a powerful tool but strongly affected by its mainly from chloroplast activity. The main mitochondrial mathematical formalization and particularly by the flux NADH producer in light conditions is the GDC system, being boundary choice (Raposo et al., 2020; Nobile et al., 2021). A consumed mainly by the malate dehydrogenase activity and not sensitivity analysis on the selected boundaries is very useful for by the electron transport chain (Bykova et al., 2014). The lack of detecting those constraints (either hard or soft constraints) that the electron transport chain in day metabolism is due to the lack generate the highest impact on the flux distribution and thus of maintenance reactions during the light period. When light is requires a thorough parameter identification. However, increased, the reactions associated with the regeneration of sensitivity analysis can also be used to detect those fluxes with RuBP are activated. This is the case of phosphoribulokinase the highest variability in relation to a given boundary, metabolic (PRK), which restores the ribulose 1,5-biphosphate. An increase model shortcomings, or fluxes that are invariable to the in the triose phosphate export activity (G3P_trans) is also boundaries. A local sensitivity analysis is presented, with a detected responding to the increase in the carbon fixation. The focus on the flux ratios defined in Table 3, which are either mitochondrial malate dehydrogenase (MDHm) is also defined empirically or based on literature and may be prone to augmented in the first place to respond to the increased uncertainty. The sensitivity analysis has also been explored for anabolic demands and in the second place given the increased the light irradiance. Details on the ranges of flux ratios explored flux through photorespiration that triggers the GCS and the for this analysis are presented in Table 6. The sensitivity analysis consequent increased consumption of NADH through MDHm has been based on the variability of the fluxes along the range of to keep the redox balance in mitochondria (Schertl and Braun, flux ratios allowed. Such variability is expressed as the 2014). The sensitivity analysis of the day metabolism also normalized summation of the slopes in each of the steps of the highlights that, when the night respiration flux in relation to range of flux ratios explored: the daily carbon fixation is increased relative to the nominal value of 0.25, the associated photosynthetic reactions (PRK and ni+1,j,z −ni,j,z i=6 Drz G3P_trans) are reduced. FVj,z = o i=1 Dvj,z (42) Regarding night metabolism (Figure 7B), the variability of the represented reactions is higher than in day metabolism. Dvj,z = max (nj,z ) − min (nj,z ) (43) When light intensity (Iu) is increased, night mitochondrial activity (MP1_n and MP2_n) is also increased to respond to In (42), the i index refers to the discrete steps in which the the maintenance reaction demands (Frantz and Bugbee, 2005). range of flux ratios has been divided (a total of six steps for each The activity of glyceraldehyde 3-phosphate dehydrogenase flux ratio) whereas index j and z refer to the reaction and the flux (GAPDc_n) responds to the increased flux towards glycolytic ratio, respectively. The parameter Drz is the increment in the flux pathways to process the carbon compounds converted during ratio expressed as a fraction and equivalent to 0.1 for all cases. night metabolic activities (Schneider et al., 2018; Gaude et al., The denominator Dvj,z is used to normalize the summation of 2018). In the dark phase, with the cyclic TCA cycle re-stored, the slopes, and it represents the range of values for a given reaction j. excess TCA intermediates and soluble sugars not used to fuel The sensitivity analysis has been carried for light and night dark metabolism are stored and reused in light metabolism for TABLE 6 Range of values explored in sensitivity analysis. Flux ratios Nominal ratio Range Description ExO2[d]: ExCO2[d] 1.22:1 [-30%, +30%] Photosynthesis rate ATPmaintenance[n]:NADPHmaintenance[n] 3:1 [-30%, +30%] Maintenance ATPmaintenance[d] :NADPHmaintenance[d] 3:1 [-30%, +30%] Maintenance ExCO2[n]:ExCO2[d] 0.25:1 [-30%, +30%] Respiration Iu 400 mmol m-2 s-1 [-30%, +30%] Irradiance The percentage is applied to the nominal ratio (i.e., for the photosynthesis rate ratio, the range explored is from 0.8:1 to 1.59:1). The whole range is split into six points. Frontiers in Plant Science 15 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 A B 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. nitrogen assimilation and biosynthesis reactions (Popova and Pinheiro De Carvalho, 1998; Igamberdiev and Eprintsev, 2016). 2us (Oin 2 − xs ) + vExO2 = 0 (44) In this way, the pool of organic carbon compounds and reducing equivalents is managed by the plants as a response to variations x^ (i + 1jk) = x^ (ijk) in the day carbon fixation efficiency and night respiration    T T activities. GAPDc_n is also positively affected when the ratio + 2U(ik) Oin 2 −x ^ (ijk) s + vExO2 s (45) V V of respiratory produced over photosynthetically consumed CO2 is increased (ExCO2[n]:ExCO2[d]) with the availability of citrate • At 10 h of simulation time, a perturbation in the process and the usage of other TCA intermediates in light photoperiods output has been introduced by adding a sudden decrease in reduced given the night metabolic increased activities as the oxygen concentration from the measurement to 20.8%. suggested by the citrate import flux. Therefore, an increased • At 16 h of simulation time, a change in the oxygen night CO2 release is corresponded by an increased glycolytic setpoint from 21% to 21.2% has been included. activity represented by GAPDC_n, also impacting the TCA intermediates and soluble sugar accumulation at night. 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 3.3 Testing the integration of multilevel offset of the O2 measurement in relation to the reference model approach to advanced control (Figure 8A). Similarly, when an abrupt perturbation is added architectures in dynamic simulations at 10 h of simulation, the gas flow is stopped to restore the oxygen concentration rapidly. Finally, when the reference is The control strategy presented based on a new approach modified from 21% to 21.2%, the gas flow is also reduced in combining the prediction capacity of MPC and the integration of order to accumulate oxygen in the growing chamber and reduce constraint-based metabolic modelling has been tested in a 24-h the tracking error (Figure 8B). The control of CO2 is achieved by dynamic simulation. To do that, different perturbations have externally injection of pure CO2 (Figure 8C). been included at different points of the simulation to test the The controller is thus demonstrated to be resilient and resiliency of the proposed control approach its capacity to smooth to overcome any of the perturbations applied as well overcome common control challenges and to analyse as on the nominal operation. This performance is due to the its versatility: reliable metabolic-based model introduced in the MPC and largely due to the integration of a disturbance specially to • At 5 h of simulation time, a plant-model mismatch has achieve offset-free control. In Figure 9, the comparison of the been introduced by adding a perturbation in the form of controller performance with and without disturbance a multiplication factor to the gas flows in (31): integration is represented. It can be observed that under Frontiers in Plant Science 16 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 A B C 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. different perturbations, when disturbance is integrated in the and surrogated models to complex genome-scale internal model of both the SSTO and the MPC, offset is reduced metabolic models. 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- 4 Conclusions model mismatches are bigger, that is, when models are less reliable than the one presented in this study. In the case of higher Modelling higher plants has major challenges to deal with. model disagreement, disturbance integration approaches should These challenges include huge metabolic changes associated be considered like integrating the error in the multistep model with the light and night photoperiods, substrate partitioning prediction (Tian et al., 2007) or integrating a moving horizon given the heterogenic requirements of the different tissues estimation (MHE) and MPC to estimate uncertainty parameters present in higher plants, organelle coordination, complex and to include them in the MPC algorithm (Huang et al., 2010). morphologies that condition the interaction with the Disturbance rejection especially in the situation of plant-model environment, and many other phenomena still not fully mismatch has been highly analysed within MPC development, understood. In this study, a multilevel model has been and certainly, the future metabolic-based controllers will need to designed with the main mechanistic phenomena that drive deal with a wide diversity of model typologies from simplified crop growth distributed into different levels in decreasing 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. Frontiers in Plant Science 17 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 order of scale length. The output of the mechanistic multilevel metabolic model development. All authors contributed to the model has been connected to a constraint-based metabolic article and approved the submitted version. model providing information of high interest about cell metabolism. The presented multilevel model offers the advantage to merge all the available information related to Funding plant growth in a structured way ensuring the solution found to be feasible for all the phenomena described in the different MELiSSA is an international consortium of 15 partners led by layers in the model hierarchy, from light reception and the European Space Agency. Its activities are governed by a biochemical conversion down to specific metabolic pathways, a Memorandum of Understanding (ESA 4000100293/10/NL/PA). feasibility that cannot be granted in only mechanistic or The MELiSSA Pilot Plant is funded from ESA contributions from metabolic-based models using stand-alone modelling Spain (main contributor), Belgium, France, Italy, and Norway, strategies. This method has been validated with experimental under Frame Contract C4000109802/13/NL/CP. Cofunding from data, integrated, and tested in a novel advanced control strategy Ministerio de Ciencia, Innovació n y Universidades, Generalitat de with promising results. Computational capabilities are no longer Catalunya, and Universitat Autònoma de Barcelona is also a major constraint, making it possible to contemplate the design acknowledged. Carles Ciurans is a PhD fellow from the POMP of control strategies that can integrate more information about Program (MELiSSA Foundation) of ESA. This work was supported the system under operation than currently done. In this study, by VILLUM FONDEN under the VILLUM Investigator Grant (no. the focus has been placed in the use of metabolic information 25920): Center for Research on Microgrids (CROM). Igor Marin de embedded in a model-predictive control providing promising Mas and Ivan Martı́nez-Mongue would like to thank The Novo results in a dynamic simulation of a growing crop chamber with Nordisk Foundation (NNF Grant numbers: NNF10CC1016517 a range of applications going from agriculture to life support and NNF14OC0009473). 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 Acknowledgments metabolic information, emerges as an opportunity to design model-based predictive controller techniques. These advanced The authors wish to thank the partners of the MELiSSA control strategies should be further explored with different levels Consortium for carefully reviewing this publication. 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 Conflict of interest increasingly available metabolic information. Author IMdM and IM are employed by Novo Nordisk Foundation Center for Sustainability. Data availability statement The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that The datasets presented in this study can be found in online could be construed as a potential conflict of interest. repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material. Publisher’s note Author contributions All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated CC is the main writer of the manuscript as well as the one organizations, or those of the publisher, the editors and the responsible for script development, simulation, and main reviewers. Any product that may be evaluated in this article, or analysis and interpretation of the results. IM and IMdM are claim that may be made by its manufacturer, is not guaranteed the main contributors to the constraint-based metabolic or endorsed by the publisher. 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 Supplementary material metabolic networks. CG and FG have contributed to the design of the mechanistic multilevel model, specially concerning the The Supplementary Material for this article can be found physical and enzymatic equations, and they have supervised the online at: https://www.frontiersin.org/articles/10.3389/ overall text. CG is the responsible of the original lettuce fpls.2022.970410/full#supplementary-material Frontiers in Plant Science 18 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 References Amitrano, C., Chirico, G. B., Pascale, S., Rouphael, Y., and Micco, V. (2020). Huang, R., Biegler, L. T., and Patwardhan, S. C. (2010). Fast offset-free nonlinear Crop management in controlled environment agriculture (CEA) systems using model predictive control based on moving horizon estimation. Ind. Eng. Chem. Res. predictive mathematical models. Sens. (Switzerland) 20 (11), 1-19. doi: 10.3390/ 49 (17), 7882–7890. doi: 10.1021/ie901945y s20113110 Igamberdiev, A. U., and Eprintsev, A. T. (2016). Organic acids: The pools of Arnold, A., and Nikoloski, Z. (2011). A quantitative comparison of Calvin- fixed carbon involved in redox regulation and energy balance in higher plants. Benson cycle models. Trends Plant Sci. 16 (12), 676–683. doi: 10.1016/ Front. Plant Sci. 7. doi: 10.3389/fpls.2016.01042 j.tplants.2011.09.004 Lee, C. P., Elsässer, M., Fuchs, P., Fenske, R., Schwarzländer, M., and Millar, A. Becker, S. A., Feist, A. M., Mo, M. L., Hannum, G., Palsson, B., and Herrgard, M. H. (2021). The versatility of plant organic acid metabolism in leaves is underpinned J. (2007). Quantitative prediction of cellular metabolism with constraint-based by mitochondrial malate-citrate exchange. Plant Cell 33 (12), 3700–3720. models: The COBRA toolbox. Nat. Protoc. 2 (3), 727–738. doi: 10.1038/ doi: 10.1093/plcell/koab223 nprot.2007.99 Liu, J., and van Iersel, M. W. (2021). Photosynthetic physiology of blue, green, Boscheri, G., Kacira, M., Patterson, L., Giacomelli, G., Sadler, P., and Furfaro, R. and red light: Light intensity effects and underlying mechanisms. Front. Plant Sci. (2012). Modified energy cascade model adapted for a multicrop lunar greenhouse 12. doi: 10.3389/fpls.2021.619987 prototype. Adv. Space Res. 50 (7), 941–951. doi: 10.1016/j.asr.2012.05.025 Louarn, G., and Song, Y. (2020). Two decades of functional-structural plant Bykova, N. V., Møller, I. M., Gardeström, P., and Igamberdiev, A. U. (2014). The modelling: Now addressing fundamental questions in systems biology and function of glycine decarboxylase complex is optimized to maintain high predictive ecology. Ann. Bot. 126 (4), 501–509. doi: 10.1093/aob/mcaa143 photorespiratory flux via buffering of its reaction products. Mitochondrion 19, Marchetti, A. G., Ferramosca, A., and Gonzá lez, A. H. (2014). Steady-state target 357–364. doi: 10.1016/j.mito.2014.01.001 optimization designs for integrating real-time optimization and model predictive Chen, X., Wang, L.c., Li, T., Yang, Q.c., and Guo, W.z. (2019). Sugar control. J. Process Control 24 (1), 129–145. doi: 10.1016/j.jprocont.2013.11.004 accumulation and growth of lettuce exposed to different lighting modes of red Marshall-Colon, A., Long, S. P., Allen, D. K., Allen, G., Beard, D. A., Benes, B., and blue LED light. Sci. Rep. 9 (1), 1–10. doi: 10.1038/s41598-019-43498-8 et al. (2017). Crops in silico: Generating virtual crops using an integrative and Cheung, C. Y. M., Williams, T. C. R., Poolman, M. G., Fell, D. A., Ratcliffe, R. G., multi-scale modeling platform. Front. Plant Sci. 8. doi: 10.3389/fpls.2017.00786 and Sweetlove, L. J. (2013). A method for accounting for maintenance costs in flux McClain, A. M., and Sharkey, T. D. (2019). Triose phosphate utilization and balance analysis improves the prediction of plant cell metabolic phenotypes under beyond: From photosynthesis to end product synthesis. J. Exp. Bot. 70 (6), 1755– stress conditions. Plant J. 75 (6), 1050–1061. doi: 10.1111/tpj.12252 1766. doi: 10.1093/jxb/erz058 Cheung, C. Y M., Poolman, M. G., Fell, D. A., George Ratcliffe, R., and Mekhalfi, M., Puppo, C., Avilan, L., Lebrun, R., Mansuelle, P., Maberly, S. C., Sweetlove, L. J. (2014). A diel flux balance model captures interactions between et al. (2014). Glyceraldehyde-3-Phosphate dehydrogenase is regulated by light and dark metabolism during day-night cycles in C3 and crassulacean acid ferredoxin-NADP reductase in the diatom asterionella Formosa. New Phytol. 203 metabolism leaves. Plant Physiol. 165 (2), 917–929. doi: 10.1104/pp.113.234468 (2), 414–423. doi: 10.1111/nph.12820 Ciurans, C., Bazmohammadi, N., Vasquez, J. C., Dussap, G., Guerrero, J. M., and Michelet, L., Zaffagnini, M., Morisse, S., Sparla, F., Pé rez-Pé rez, M. E., Francia, Godia, F. (2021). Hierarchical control of space closed ecosystems – expanding F., et al. (2013). Redox regulation of the Calvin-Benson cycle: Something old, microgrid concepts to bioastronautics. IEEE Ind. Electron. Mag. 15 (2), 2–13. something new. Front. Plant Sci. 4. doi: 10.3389/fpls.2013.00470 doi: 10.1109/MIE.2020.3026828 Miginiac-Maslow, M., and Lancelin, J. M. (2002). Intrasteric inhibition in redox Corpas, F. J., and Barroso, J. B. (2014). NADPH-generating dehydrogenases: signalling: Light activation of NADP-malate dehydrogenase. Photosyn. Res 72 (1), Their role in the mechanism of protection against nitro-oxidative stress induced by 1-12. doi: 10.1023/A:1016099228450 adverse environmental conditions. Front. Environ. Sci. 2. doi: 10.3389/ fenvs.2014.00055 Morcuende, R., Krapp, A., Hurry, V., and Stitt, M. (1998). Sucrose-feeding leads to increased rates of nitrate assimilation, increased rates of a-oxoglutarate Dong, C., Fu, Y., Xie, B., Wang, M., and Liu, H. (2017). Element cycling and synthesis, and increased synthesis of a wide spectrum of amino acids in tobacco energy flux responses in ecosystem simulations conducted at the Chinese lunar leaves. Planta 206, 394–409. doi: 10.1007/s004250050415 palace-1. Astrobiology 17 (1), 78–86. doi: 10.1089/ast.2016.1466 Morgan, J. A., and Rhodes, D. (2002). Mathematical modeling of plant metabolic Ebrahim, A., Lerman, J. A., Palsson, B. O., and Hyduke, D. R. (2013). COBRApy: pathways. Metab. Eng. 4 (1), 80–89. doi: 10.1006/mben.2001.0211 COnstraints-based reconstruction and analysis for Python. BMC Syst. Biol. 7, 1-6. doi: 10.1186/1752-0509-7-74 Nikolov, N. T., Massman, W. J., and Schoettle, A. W. (1995). Coupling biochemical and biophysical processes at the leaf level: An equilibrium Eckart, P. (1995). Spaceflight life support and biospherics. (Dordrecht: Springer), photosynthesis model for leaves of C3 plants. Ecol. Model. 80 (2–3), 205–235. 1–444. doi: 10.1016/0304-3800(94)00072-P Farquhar, G. D., Caemmerer, S.v., and Berry, J. A. (1980). A biochemical model Nobile, M. S., Coelho, V., Pescini, D., and Damiani, C. (2021). Accelerated global of photosynthetic CO2 assimilation in leaves of C3 species. Planta 149 (1), 78–90. sensitivity analysis of genome-wide constraint-based metabolic models. BMC doi: 10.1007/BF00386231 Bioinf. 22, 1-17. doi: 10.1186/s12859-021-04002-0 Frantz, J. M., and Bugbee, B. (2005). Acclimation of plant populations to shade: Ocampo-Martinez, C., Barcelli, D., Puig, V., and Bemporad, A. (2012). Photosynthesis, respiration, and carbon use efficiency. J. Am. Soc. Hortic. Sci. 130 Hierarchical and decentralised model predictive control of drinking water (6), 918–927. doi: 10.21273/jashs.130.6.918 networks: Application to Barcelona case study. IET Control Theory Appl. 6 (1), Fulget, N., Poughon, L., Richalet, J., and Lasseur, C. (1999). MELISSA : Global 62–71. doi: 10.1049/iet-cta.2010.0737 control strategy of the artificial ecosystem by using first principles models of the Pannocchia, G., and Bemporad, A. (2007). Combined design of disturbance compartments. Adv. Space Res. 24 (3), 397–405. doi: 10.1016/S0273-1177(99) model and observer for offset-free model predictive control. IEEE Trans. Automat. 00490-1 Control 52 (5), 1109–1113. doi: 10.1109/TAC.2007.899096 Gakière, B., Hao, J., Bont, L., Pé triacq, P., Nunes-Nesi, A., and Fernie, A. R. Peng, B., Guan, K., Tang, J., Ainsworth, E. A., Asseng, S., Bernacchi, C. J., et al. (2018). NAD + biosynthesis and signaling in plants. Crit. Rev. Plant Sci. 37 (4), (2020). Towards a multiscale crop modelling framework for climate change 259–307. doi: 10.1080/07352689.2018.1505591 adaptation assessment. Nat. Plants 6 (4), 338–348. doi: 10.1038/s41477-020- Gaude, E., Schmidt, C., Gammage, P. A., Dugourd, A., Blacker, T., Chew, S. P., 0625-3 et al. (2018). NADH shuttling couples cytosolic reductive carboxylation of Popova, T. N., and Pinheiro De Carvalho, M. A. A. (1998). Citrate and isocitrate glutamine with glycolysis in cells with mitochondrial dysfunction. Mol. Cell 69 in plant metabolism. Biochim. Biophys. Acta 1364, 307–325. doi: 10.1016/S0005- (4), 581–593.e7. doi: 10.1016/j.molcel.2018.01.034 2728(98)00008-5 Gitelson, J., and Lisovsky, G. (2002). Man-made closed ecological systems Poughon, L., Farges, B., Dussap, C. G., Godia, F., and Lasseur, C. (2009). (London: CRC Press). Simulation of the MELiSSA closed loop system as a tool to define its integration Gomes, C., Dal’molin, O., Quek, L.-E., Saa, P. A., Nielsen, L. K., Nikoloski, Z., strategy. Adv. Space Res. 44 (12), 1392–1403. doi: 10.1016/j.asr.2009.07.021 et al. (2015). A multi-tissue genome-scale metabolic modeling framework for the Poulet, L. (2018). Developpement de modé les physiques pour comprendre la analysis of whole plant systems. Front. Plant Sci 6. doi: 10.3389/fpls.2015.00004 croissance des plantes en environment de gravité ré duite pour des applications dans Harley, P. C., Thomas, R. B., Reynolads, J. F., and Strain, B. R. (1992). Modelling Les systé mes suuport-vie (Clermont-Ferrand: Université Clermont-Auvergne). photosynthesis of cotton grown in elevated CO2. Plant Cell Environ. 15 (3), 271– Poulet, L., Dussap, C. G., and Fontaine, J. P. (2020). Development of a 282. doi: 10.1111/j.1365-3040.1992.tb00974.x mechanistic model of leaf surface gas exchange coupling mass and energy Frontiers in Plant Science 19 frontiersin.org Ciurans et al. 10.3389/fpls.2022.970410 balances for life-support systems applications. Acta Astronaut. 175, 517–530. Taniguchi, M., and Miyake, H. (2012). Redox-shuttling between chloroplast and doi: 10.1016/j.actaastro.2020.03.048 cytosol: Integration of intra-chloroplast and extra-chloroplast metabolism. Curr. Poulet, L., Fontaine, J. P., and Dussap, C. G. (2018). A physical modeling Opin. Plant Biol. 15 (3), 252–260. doi: 10.1016/j.pbi.2012.01.014 approach for higher plant growth in reduced gravity environments. Astrobiology 18 Tatjewski, P. (2017). Offset-free nonlinear model predictive control with state- (9), 1093–1100. doi: 10.1089/ast.2017.1804 space process models. Arch. Control Sci. 27 (4), 595–615. doi: 10.1515/acsc-2017-0035 M. Raposo, P. Ribeiro, S. Sé rio, A. Staiano and A. Ciaramella (Eds.) (2020). Thornley, J. H. M. (2002). Instantaneous canopy photosynthesis: Analytical Computational intelligence methods for bioinformatics and biostatistics (Cham: expressions for sun and shade leaves based on exponential light decay down the Springer International Publishing). Available at: https://link.springer.com/10.1007/ canopy and an acclimated non-rectangular hyperbola for leaf photosynthesis. Ann. 978-3-030-34585-3. Bot. 89 (4), 451–458. doi: 10.1093/aob/mcf071 Sasidharan, S. (2012). Modélisation de la croissance des plantes supérieures pour Thornley, J., and Johnson, I. (1980). Plant and crop modelling: A mathematical les systèmes de support-vie : modèle métabolique de la feuille de laitue considérant approach to plant and crop phisiology (Oxford: The Blackburn Press). la conversion d’énergie et le métabolisme central du carbone [PhD dissertation]. Tian, X., Wang, P., Huang, D., and Chen, S. (2007). Offset-free multistep [Clermont-Ferrand]: Université Blaise Pascal. nonlinear model predictive control under plant-model mismatch. Int. J. Adapt. Scattolini, R. (2009). Architectures for distributed and hierarchical model Control Signal Process. 21, 731–744. doi: 10.1002/acs predictive control - a review. J. Process Control 19 (5), 723–731. doi: 10.1016/ Vasquez, J., Guerrero, J., Miret, J., Castilla, M., and Vicuña, L. (2010). j.jprocont.2009.02.003 Hierarchical control of intelligent microgrids. IEEE Ind. Electron. Mag. 4 (4), 23– Schertl, P., and Braun, H. P. (2014). Respiratory electron transfer pathways in 29. doi: 10.1109/MIE.2010.938720 plant mitochondria. Front. Plant Sci. 5. doi: 10.3389/fpls.2014.00163 Xiao, Y., Chang, T., Song, Q., Wang, S., Tholen, D., Wang, Y., et al. (2017). Schneider, M., Knuesting, J., Birkholz, O., Heinisch, J. J., and Scheibe, R. (2018). EPlant for quantitative and predictive plant science research in the big data era–lay Cytosolic GAPDH as a redox-dependent regulator of energy metabolism. BMC the foundation for the future model guided crop breeding, engineering and Plant Biol. 18 (1), 1–14. doi: 10.1186/s12870-018-1390-6 agronomy. Quant. Biol. 5 (3), 260–271. doi: 10.1007/s40484-017-0110-9 Shameer, S., Ratcliffe, R. G., and Sweetlove, L. J. (2019). Leaf energy balance Yu, C., Liang, Z., and Huang, A. H. C. (1984). Glyoxylate transamination in requires mitochondrial respiration and export of chloroplast NADPH in the light. intact leaf peroxisomes. Plant Physiol. 75 (1), 7–12. doi: 10.1104/pp.75.1.7 Plant Physiol. 180 (4), 1947–1961. doi: 10.1104/pp.19.00624 Zhou, W., Liang, X., Li, K., Dai, P., Li, J., Liang, B., et al. (2021). Metabolomics Soussana, J.-F. X., Graux, A.-I., and Tubiello, F. N. (2010). Improving the use of analysis reveals potential mechanisms of phenolic accumulation in lettuce (Lactuca modelling for projections of climate change impacts on crops and pastures. J. Exp. sativa l.) induced by low nitrogen supply. Plant Physiol. Biochem. 158, 446–453. Bot. 61 (8), 2217–2228. doi: 10.1093/jxb/erq100 doi: 10.1016/j.plaphy.2020.11.027 Tan, X. L. J., and Cheung, C. Y. M. (2020). A multiphase flux balance model Zhou, J., Wang, J. Z., Hang, T., and Li, P. P. (2020). Photosynthetic characteristics and reveals flexibility of central carbon metabolism in guard cells of C3 plants. Plant J. growth performance of lettuce (Lactuca sativa l.) under different Light/Dark cycles in 104 (6), 1648–1656. doi: 10.1111/tpj.15027 mini plant factories. Photosynthetica 58 (3), 740–747. doi: 10.32615/ps.2020.013 Frontiers in Plant Science 20 frontiersin.org