Coupled Gas-Exchange Model for C4 Leaves Comparing Stomatal Conductance Models

Plant simulation models are abstractions of plant physiological processes that are useful for investigating the responses of plants to changes in the environment. Because photosynthesis and transpiration are fundamental processes that drive plant growth and water relations, a leaf gas-exchange model that couples their interdependent relationship through stomatal control is a prerequisite for explanatory plant simulation models. Here, we present a coupled gas-exchange model for C4 leaves incorporating two widely used stomatal conductance submodels: Ball–Berry and Medlyn models. The output variables of the model includes steady-state values of CO2 assimilation rate, transpiration rate, stomatal conductance, leaf temperature, internal CO2 concentrations, and other leaf gas-exchange attributes in response to light, temperature, CO2, humidity, leaf nitrogen, and leaf water status. We test the model behavior and sensitivity, and discuss its applications and limitations. The model was implemented in Julia programming language using a novel modeling framework. Our testing and analyses indicate that the model behavior is reasonably sensitive and reliable in a wide range of environmental conditions. The behavior of the two model variants differing in stomatal conductance submodels deviated substantially from each other in low humidity conditions. The model was capable of replicating the behavior of transgenic C4 leaves under moderate temperatures as found in the literature. The coupled model, however, underestimated stomatal conductance in very high temperatures. This is likely an inherent limitation of the coupling approaches using Ball–Berry type models in which photosynthesis and stomatal conductance are recursively linked as an input of the other.


Introduction
Leaf gas-exchange includes the processes of CO 2 assimilation and water vapor exchange by plant leaves. It is one of the most important processes for life on Earth as it provides carbohydrate for food and oxygen for respiration practically for all organisms. Because plant growth depends on photosynthesis, it is an essential building block of plant simulation models. Photosynthesis models range in complexity from correlative models based on radiation use efficiency where carbon assimilation is proportional to total irradiance absorbed by leaf surfaces, to models based on enzyme kinetics [1]. Integration of more mechanistic photosynthesis model has been a critical aspect in crop modeling to better understand and predict crop productivity under dynamic environments [2,3].
A coupled approach to modeling photosynthesis, stomatal conductance, and transpiration simultaneously for C 3 plants has been presented by a number of studies [4][5][6][7][8]. This approach usually combines the FvCB (Farquhar-von Caemmerer-Berry) C 3 photosynthesis model [9] with a model of

Gas-Exchange Model Calibration
The performance of calibrated parameters was evaluated by Willmott's refined index of agreement (d r ) and Nash-Sutcliffe model efficiency coefficient (NSE) for the two variants of gas-exchange model ( Figure 1). d r of the model with Ball-Berry stomatal conductance submodel (BB) was 0.879 for net photosynthesis rate (A n ) and 0.804 for stomatal conductance (g s ). d r of the other model with Medlyn submodel (MED) was 0.881 for net photosynthesis rate (A n ) and 0.820 for stomatal conductance (g s ). NSE of the BB model was 0.941 for A n and 0.798 for g s . NSE of the MED model was 0.937 for A n and 0.796 for g s . d r values were close to 1, meaning our models were calibrated adequately for further analysis carried in the following sections, especially when it was achieved by calibrating only four parameters, while most other parameters came from existing models or literature (Table A2). The two variants of the gas-exchange model with different selection of the stomatal conductance submodel did not show a clear difference in their performance within the range of input given by our experimental dataset.
Although the difference between two gas-exchange models only came from how stomatal conductance was calculated in submodel, the pooled nature of our calibration process produced a slightly different set of parameters used by the other parts of model, such as nitrogen submodel shared between the two variants. For example, with baseline leaf nitrogen content (N 0 ) and steepness of nitrogen response curve (s) parameters, BB got 0.371 and 4.470, while MED got 0.315 and 3.912, respectively. To ensure their difference did not hamper comparison between two stomatal conductance models, we did a further sensitivity analysis on the two other parameters in Appendix C and confirmed the difference was negligible in regards with A n response ( Figure A7).

Stomatal Conductance Model Comparison
g s predicted by MED was generally higher than g s from BB ( Figure 2). g s from both models were in a similar range when relative humidity (RH) was around 65 % to 90 %, while maximum g s was much higher in MED when RH was saturated. When RH went lower than 50 %, gs from BB dropped rapidly and mostly converged to the lower bound (g 0 BB ), almost shutting down transpiration. In MED, decrease of g s along RH gradient was more gradual, and its value usually remained higher than the lower bound (g 0 MED ) to ensure a certain amount of transpiration keeps occurring. Stomatal conductance (g s ) estimated by two stomatal conductance models over a range of atmospheric CO 2 concentration (C a ) from 10 µbar to 1500 µbar at multiple levels of relative humidity (RH). The graph was plotted against resultant intercellular CO 2 concentration (C i ). Air temperature (T a ) was 32 • C and irradiance (I) was 2000 µmol quanta m −2 s −1 .
As a result, A n from the gas-exchange model using BB decreased more rapidly, while the counterpart using MED showed a much gentle response as RH went down at the same level of concentration of atmospheric CO 2 (C a ). The rate of decrease in A n was 75 % with BB and 20 % with MED when RH was dropped from 80 % to 20 %. However, there was little difference in the curvature of A n response between the two models in terms of intercellular CO 2 concentration (C i ), suggesting the difference was mostly due to a change in supply function of the curve ( Figure 3).  . Net photosynthesis rate (A n ) estimated by two stomatal conductance models over a range of atmospheric CO 2 concentration (C a ) from 10 µbar to 1500 µbar at two levels of relative humidity (RH). The graph was plotted against resultant intercellular CO 2 concentration (C i ). Air temperature (T a ) was 32 • C and irradiance (I) was 2000 µmol quanta m −2 s −1 .
A similar response can be observed when A n was plotted against a range of air temperature (T a ) from 0 • C to 50 • C (Figure 4). At lower RH, MED was able to maintain g s to a certain level and therefore keep A n from collapsing. The shifting of optimal temperature towards a higher regime with lower RH came from increased cooling effect by higher water loss in drier condition.  In turn, leaf temperature (T l ), which is adjusted by energy balance equation involving latent heat flux mainly driven by leaf transpiration, showed a clear difference between the two models. With BB, lower RH had stomata almost closed down and thus not able to cool down leaf temperature with latent cooling. On the contrary, MED implied higher latent cooling under lower RH due to stronger gradient of water vapor pressure formed between air and inside the leaf.

Leaf Nitrogen Deficiency
A n generally decreased as leaf nitrogen content (N) reduced and the rate of decrease accelerated when nitrogen was more limited as represented by logarithmic curves. Yet, the strength of response was not constant and varied depending on environmental conditions. A n was more decreased with lower RH, but the difference diminished when N went below 0.5 g m −2 or relative leaf nitrogen content (N p ) was less than 1 % assuming specific leaf area (SLA) was 200 cm 2 g −1 . The rate of A n decrease did not change much with high atmospheric CO 2 (C a ) and only had more negative effect when C a was below 400 µmol mol −1 . Response to T a under nitrogen stress was nonlinear that the decrease was more significant when T a was moving away from optimal temperature. The optimal temperature, where a peak A n could be achieved, slightly increased with more N available; thus, the maximum A n itself also increased due to more favorable biochemical reactions with higher temperature and N (Figure 5c,g). The slope of A n decrease by N stress was steeper under higher irradiance (I) and the difference between the levels of I gradually diminished as N approaching a minimum (Figure 5a,e). Overall, there was no clear difference in terms of nitrogen response between BB and MED models.

Leaf Water Status
A n decreased as bulk leaf water potential (Ψ v ) reduced, but the rate of decrease did not monotonically change as in the case of nitrogen stress. Generally under a greater water stress with lower water potential, stress response represented by A n reduction tapered off and formed a logistic response curve. A n was more decreased with lower RH, but the difference diminished when Ψ v kept decreasing. In addition, note that with higher RH, leaf was able to sustain maximum A n even under mild water deficit which would have led to a noticeable reduction in A n under lower RH. For example, at −0.4 MPa, A n under 80 % of RH did not decrease much, whereas A n under 40 % of RH saw almost 60 % reduction with BB ( Figure 5b) and 20 % reduction with MED ( Figure 5f). The decrease of A n was more consistent with all range of C a compared to nitrogen stress response. For −1.5 MPa and below, the rate of A n decrease became almost identical regardless of C a . With higher C a , maximum A n was sustained for larger range of Ψ v similar to the RH response. Water response to T a was also nonlinear, but the difference wore off with lower Ψ v and no single optimal temperature was clear to be found (Figure 5d,h). Overall response to T a under water stress was similar between BB and MED, except BB exhibited much higher optimal temperature than MED in a mild-to-severe stress level, indicated by −1.0 MPa. Overall, BB was more sensitive to RH changes under water stressed conditions than MED (Figure 5d,f).

Response to Relative Humidity
At high RH, A n remained relatively stable within its optimal range when both N and Ψ v kept high for less stress (Figure 6a,g). An area of this region located in the upper right side of contour plot shrunk and so did the range of non-limiting N and Ψ v as RH went down. In other words, with lower RH, A n became more sensitive to both nitrogen and water stress factors. The reduction of A n for lower RH was extremely strong with BB ( Figure 6d) compared to MED (Figure 6j). No other variables showed such drastic difference between BB and MED in the comparison.
Note that, with low N and high Ψ v , the sensitivity of A n was mostly governed by the change of N, thus limited by nitrogen. On the other hand, with high N and low Ψ v , the most of sensitivity came from the change of Ψ v , thus limited by water. Then, with both low N and Ψ v , A n was also largely driven by Ψ v unless N dropped down to a very low range.

Response to Atmospheric CO 2
The upper-right region of non-limiting A n existed with a range of C a above 400 µbar. This region vertically expands further down to cover lower Ψ v under higher C a , as shown by comparing ambient CO 2 (Figure 6e,k) with elevated CO 2 (Figure 6b,h), indicating possible alleviation of water stress by elevated CO 2 concentration. A n remained relatively stable until Ψ v reached down to −1.0 MPa under 800 µbar of C a , whereas A n started decreasing faster only after −0.5 MPa under 400 µbar of C a . Below these boundaries, A n was mostly limited by water only.

Response to Irradiance
The region of non-limiting A n remained relatively stable with a range of I (Figure 6c,i) and even expanded further down when I was lower, although the magnitude of A n was much smaller (Figure 6f,l). In other words, under shaded condition with less light, nitrogen and water became less relevant because A n had to be much smaller. Below this region, A n was again mostly limited by water only unless N deficiency was extremely strong. Ball

Performance of Stomatal Conductance Models
We did not find much difference in performance between two stomatal conductance submodels BB and MED when used for fitting calibration dataset ( Figure 1). It is known that two models have equal predictive strength for non-extreme environmental conditions where the data are usually collected for calibration and validation [25]. Results started deviating from each other under more extreme conditions, such as low C a and low RH, which are also close to the range where gas-exchange instruments often find difficulties in accurate measurements (Figures 2-5b,f). Under low RH, especially below 50 %, g s BB had much higher rate of decrease and almost converged to g 0 BB in the end, whereas g s MED maintained higher than g 0 MED with smooth and gentle transitions most of the time ( Figure 2). Seemingly degenerate behavior of g s BB might be an overlooked effect of g 0 BB as a parameter estimated from regression rather than a directly measured value [26]. A clearly different response from each model warrants caution, especially when applying models to extreme conditions, such as concurrent exposure to high temperature and high vapor pressure deficit (VPD), likely encountered in future climate projections [27].

Stress Response to Elevated CO 2
Nitrogen stress consistently posed a negative impact on A n although the degree of reduction may vary depending on the severity of stress and corresponding environmental conditions (Figure 5a,c,e,g). While water stress also reduced A n in most of the time (Figure 5b,d,f,h), we observed in some conditions that the impact of water stress greatly diminished, for instance, when RH was higher than 60 % (Figure 6a,g) or C a was higher than 400 µmol mol −1 (Figure 6b,h). There is increasing evidence that elevated CO 2 concentration alleviates water stress [28][29][30][31][32][33][34]. Our simulation result confirms this positive effect of elevated CO 2 by showing high A n sustained in a wider range of Ψ v under high C a due to reduced g s preserving water loss via transpiration. The effect was mostly pronounced under a mild water stressed condition since maximum A n did not change much in the absence of water stress [35]. Alleviation strength was dependent on leaf nitrogen supply such that even very mild water stress could not be overcome by high C a under low N [33].

Interactions between Nitrogen Deficiency and Water Stress
When both nitrogen and water stress were imposed simultaneously, the net effect of stress may vary depending on the relative degree of each stress factor and environmental variables. In the contour plots of interactive stress effects, we can identity three types of contours ( Figure 5). Mostly vertical lines with very steep slopes indicate A n could be only improved by moving across a horizontal gradient, conferring a dominant N sensitivity. Likewise, mostly horizontal lines with a slope close to zero indicate A n could be only improved by moving across a vertical gradient, conferring a dominant Ψ v sensitivity. Other curves from round contour lines indicate A n are affected by both stress factors. With that in mind, looking back stress interaction result figures gives an insight that a vast area of N and Ψ v grid are covered by horizontal lines; thus, they are mostly under water limiting conditions. Nitrogen limiting conditions are, on the other hand, exhibited within an area where N is low or Ψ v is high. Such distinctions could be also found in literature where the negative impact of nitrogen stress could be inflated by an existence of drought stress [30,[36][37][38].

Response of Transgenic C 4 Leaves to Temperature
As an extended exercise for testing model response against real conditions, we applied our model to simulate a wide range of temperature response of transgenic C 4 plants, Flaveria bidentis, with varying amount of Rubisco reported by Kubien et al. [39]. By using a set of slightly modified parameters, as described in Appendix D, we were able to replicate overall patterns of A n ( Figure 7d) and C i ( Figure 7f) over a wide range of T l using MED. The performance of BB was close to MED, except that C i , under low T a with more reduced Rubisco (AR2), showed an abrupt change due to g s that went too low, reaching the lower bound g 0 BB (Figure 7a,c).

Ball-Berry (BB)
Tl (°C)  Table A3. Other parameters remained the same as listed in Table A2. Solid curves represent model estimation for each treatment. Dots indicate data points digitized from the result of original experiment.
Simulated g s closely followed the actual measurements until around 35 • C where model and observation started showing disagreement (Figure 7b,e). In the experiment, g s kept increasing further with higher T l , whereas our model estimated g s should decrease when A n was not able to keep up under excessive heat stress. Such decreasing patterns occurred regardless of stomatal conductance submodels. This behavior is indeed contrary to observations where transpiration rate continually increases with higher leaf temperature, implying A n inhibition is not caused by stomatal closure; thus, A n should decouple from g s for some species [40,41]. However, this mechanism is not easily applicable to many if not all commonly used gas-exchange models relying on empirical relationship between A n and g s , as explained in Appendix A.5. Alternative stomatal conductance models with mechanistic understandings on stomata responses to humidity and temperature could be useful for simulating field conditions where temperature fluctuates more than humidity [42,43]. In the meantime, more attention should be paid to the current results of gas-exchange models under high temperature, which is especially critical when simulating climate change scenarios.

Model Structure
The gas-exchange model consists of a number of smaller submodels describing different aspect of biochemical and physical processes coupled to each other following the structure of Kim and Lieth [7]. The C 4 photosynthesis is coupled with a stomatal conductance model via net photosynthesis rate (A n ) and also interacts with energy balance model to establish leaf temperature (T l ) [44]. We tested two stomatal conductance models, Ball-Berry model (BB) and Medlyn model (MED), combined with the rest of submodels remained the same, yielding two variations of the gas-exchange model. More details about the model specification are given in Appendix A.

Model Calibration
Most of parameter values came from existing literatures (Table A2). Some photosynthetic parameters were obtained from Pioneer hybrid 3733 maize (Zea mays) grown under Soil-Plant-Atmosphere Research (SPAR) chambers located at Beltsville, MD, USA, in 2002 [45]. Only four parameters were specifically calibrated for modeling experiments presented in this paper. Two are nitrogen-related parameters: baseline leaf nitrogen content (N 0 ) and steepness of nitrogen response curve (s). The other two are related to stomatal conductance: lower bound of stomatal conductance (g 0 ) and sensitivity of stomatal conductance (g 1 ). For comparison of overall gas-exchange response between two stomatal conductance models, calibration was separately done for the two models: Ball-Berry (BB) and Medlyn (MED) models.
N 0 and s calibrated for BB and MED were then pooled together to provide average parameter values used by the both models. To make sure a slight divergence in the parameter value did not make an impact on overall model response, we also conducted a sensitivity analysis on N 0 and s for a range of possible values. As g 0 and g 1 should involve inherent differences between two stomatal conductance models, a separate set of parameter values were used accordingly. BB model relied on g 0 BB and g 1 BB . MED model relied on g 0 MED and g 1 MED . The optimization method used for calibration was differential evolution algorithm [46].
The experimental dataset for calibration was collected from a growth chamber experiment conducted with Pioneer hybrid 34N43 maize in 2005 at Beltsville, MD, USA, as used in Kim et al. [44,47]. Spot measurements of gas-exchange were recorded with LI-COR LI-6400XT when plants were at the onset of reproductive stage. Three nitrogen levels at 0 kg ha −1 , 50 kg ha −1 and 200 kg ha −1 in total were applied two times. For leaf nitrogen content (N), SPAD measurements were instead collected with Konica Minolta SPAD-502 and converted to N, as described in Appendix B. The fitness of model calibration was evaluated by Willmott's refined index of agreement d r as defined in Equation (1) [48] and Nash-Sutcliffe modeling efficiency coefficient (NSE) as defined in Equation (2) [49]. y i is an observed value for the variable of interest under a specific environmental condition ordered by i.ŷ i is an estimation by the model for the same input condition as y i .ȳ is the mean of observed values. Net photosynthetic rate (A n ) and stomatal conductance (g s ) were two variables selected for the evaluation.

Model Comparison
We compared two variants of the gas-exchange model depending on which submodel was used for calculating stomatal conductance: Ball-Berry (BB) and Medlyn (MED). Net photosynthetic rate (A n ) was calculated over a range of values for specific environmental input variables, atmospheric CO 2 and air temperature (T a ). C a ranged from 0 µbar to 1500 µbar and T a ranged from −10 • C to 50 • C. Other environmental variables remained constant as default values (Table 1). Since stomatal conductance can be more directly related to the relative humidity of surrounding air (RH) than other variables, a response curve of A n was separately obtained for ten levels of RH from 0 % to 100 %.

Model Response
We simulated nitrogen and water stress factors by adjusting relevant parameter values. For nitrogen stress, leaf nitrogen content (N) was changed from 0 g m −2 to 2 g m −2 to impose stress when N was lower. For water stress, bulk leaf water potential (Ψ v ) was changed from 0 MPa down to −3 MPa to impose stress when Ψ v was lower. The range of values were deliberately chosen by first selecting biophysically feasible extremes, 0 g m −2 for N and 0 MPa for Ψ v , then selecting a value for the other end where output variable exhibits a clear convergence. When simulating a response for each stress factor, the other stress factor assumed to be not limiting. The response curve of net photosynthetic rate (A n ) was obtained for multiple environmental input variables including relative humidity (RH), atmospheric CO 2 (C a ), air temperature (T a ), and irradiance (I).
In order to observe interaction effects between the two stress factors, contour plot of A n was obtained by applying the same range of two stress factors simultaneously. For each environmental input variable, four different levels of the value were chosen to show a varying contrast as the input changes. RH changed to 30, 50, 70 and 90 %. C a changed to 200, 400, 600 and 800 µbar. I changed to 500, 1000, 1500 and 2000 µmol m −2 s −1 .

Conclusions
We compared how two stomatal conductance models, Ball-Berry (BB) and Medlyn (MED), perform when coupled with the C 4 photosynthesis model of von Caemmerer [12]. Traditionally, BB has been widely used in leaf gas-exchange modeling research. More recently proposed MED approach is founded upon BB approach to provide physiological underpinnings of empirical parameters and has been considered as an alternative to BB. Our results confirmed that the performance of two models were comparable in a wide range of environmental conditions, yet could deviate substantially in low humidity conditions.
We further tested applicability of the model by replicating the behavior of transgenic C 4 leaves under moderate temperatures found in the literature. The coupled model, however, underestimated stomatal conductance in very high temperatures, presumably due to an inherent limitation of the coupling approaches using Ball-Berry type models in which photosynthesis and stomatal conductance are recursively linked as an input of the other.
We were able to reach these findings thanks to a novel modeling framework written in Julia programming language used for model development and analysis [50]. Our modeling approach is extensible and can be a useful means for studying the ecophysiology of C 4 plants, including staple and energy crops, such as maize, sorghum, and switchgrass.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: BB Ball-Berry stomatal conductance model [

Appendix A. Gas-Exchange Model Specification
Appendix A. 1

. C 4 Photosynthesis
The biochemical demand for CO 2 assimilation in C 4 leaves was adapted from an existing C 4 photosynthesis model [12]. The rate of net CO 2 assimilation (A n ) was represented by the minimum of enzyme limited (A c ) and electron transport limited (A j ) CO 2 assimilation rates with curvature factor β.
The transition between A c and A j was calculated using a hyperbolic minimum where minh{a, b, c} is equivalent to taking the lower root of quadratic equation cx 2 − (a + b)x + ab = 0 with curvature factor c interpreted as a parameter of co-limitation [51,52]. A c can be approximated by the minimum of phosphoenolpyruvate carboxylase (PEPC) activity (A c 1 ) and Rubisco activity (A c 2 ) taking into account the bundle-sheath leakage and mitochondrial respiration.
V p is the rate of C 4 carboxylation and assumed to be limited either by PEPC activity or PEP regeneration. g bs is the bundle sheath conductance to CO 2 and C m is the mesophyll CO 2 partial pressure.
V cmax is the maximum rate of Rubisco carboxylation (µmol m −2 s −1 ). R d is mitochondrial respiration in the light (µmol CO 2 m −2 s −1 ) and half of its value was assumed a mesophyll component, R m .
V pmax is the maximum PEP carboxylation rate, K p is the Michaelis-Menton constant for CO 2 of PEPC, and V pr is PEP regeneration rate. Provided that the resistance for CO 2 from intercellular spaces to mesophyll cells is negligible, C m can be estimated from CO 2 partial pressure of the air (C a ) after taking account of total leaf resistance to CO 2 (r v c ) under a given A n .
The ratio between diffusion coefficient for water vapor (D w ) and CO 2 (D c ) is used for converting from stomatal conductance (g s ) and boundary layer conductance (g b ) in terms of CO 2 into resistance in terms of water vapor (r s c , r b c ). We assume g s is subject to still air and thus raised to the power of 1, whereas g b is subject to convective air with laminar flow and thus raised to the power of 2 3 [53,54]. C m was solved numerically by using bisection method because its dependence on A n which is required in the calculation of g s forms a cyclic dependency difficult to be solved analytically.
The assimilation rate limited by electron transport (A j ) can be approximated similarly to A c by the minimum of electron transport limited rates in C 4 and C 3 cycles.
Total rate of electron transport (J) was modeled using a non-rectangular hyperbola which can be described with hyperbolic minimum (minh). x is a partitioning factor of J.
θ is curvature of response of electron transport to photosynthetically active radiation (PAR) and J max is the maximum rate of electron transport. I 2 is effective radiation absorbed by Photosystem II (PSII) [7,55].
I is the incident light in photosynthetic photon flux density (PPFD) and I a is PPFD absorbed by the leaf. α is leaf absorptance in PAR and assumed 1 − δ where δ is the proportion of the incident light scattered by the leaf surface. f is the spectral correction factor. PAR assumes the wave band of solar radiation from 400 nm to 700 nm.
The temperature dependence of V pmax , V cmax , and R d was approximated by the Arrhenius equation (k T A ) normalized at a base temperature.
E a is the activation energy varies by process, T l is the leaf temperature and T b is the base temperature assumed 25 • C. T l k and T b k are corresponding absolute temperatures in Kelvin (K). R is the universal gas constant. Each temperature dependent parameter has a value at 25 • C.
The temperature dependence of K p and V pr was assumed to be Q 10 of 2.0.
The temperature dependence of J max was modeled using a peaked function [56].
E a and R are as defined above. H is the curvature parameter determining the rate of decrease above the peak temperature and S is the entropy factor.
Nitrogen dependence of V pmax , V cmax , and J max was modeled by a logistic function for the current leaf nitrogen content (N) [20]. N 0 is a baseline value assumed for the leaf nitrogen content and s is the steepness of nitrogen response curve.
Temperature dependence (k T ) and nitrogen dependence (k N ) are multiplicative limiting factors. For example, electron transport rate (V pmax ) is scaled from V pmax 25 by calculating V pmax = V pmax 25 · k T A [T l , E ap ] · k N with the current leaf temperature (T l ) and an activation energy for PEPC (E ap ).

Appendix A.2. Boundary Layer Conductance
Boundary layer conductance to water vapor (g b ) was derived from convective heat conductance of the leaf surface (g H ). Corresponding Nusselt number (Nu) and Reynolds Number (Re) were derived by assuming forced convection of streamline flow on flat plates [57]. u is wind speed (m s −1 ) and d is characteristic dimension of the leaf determined by leaf width W (cm) [53]. D m is kinematic viscosity of the air and D h is thermal diffusivity of the air. D m and D h were assumed temperature independent at 20 • C.
g H (m s −1 ) is then scaled and converted to g h (mmol m −2 s −1 ) using the Gas Laws [57]. T a k is the absolute temperature of air (K).
g b was then derived from g h using a diffusion ratio between water vapor and heat transfer under forced convection [53,54]. Additionally, g b was normalized to the current atmospheric pressure (P a ) to make its units explicitly expressed in vapor pressure gradient rather than in unitless mole fraction.
We used two approaches for modeling stomatal conductance (g s ). g s BB is stomatal conductance following Ball-Berry (BB) model [10] and g s MED follows Medlyn (MED) model [11].
is residual stomatal conductance to water vapor at the light compensation point and g 1 BB is the empirical coefficient for the sensitivity of g s to other variables. h s is relative humidity at the leaf surface (as a fraction), A n is net photosynthesis, and P a is the partial pressure of the air.
f Ψ v is a weight factor that adjusts stomatal conductance to the bulk leaf water potential (Ψ v ) where Ψ f is a reference potential and s f is a sensitivity parameter [8].
g s MED has a slightly different form compared to g s BB that uses vapor pressure deficit at the leaf surface (D s ) instead of h s . g 0 MED (mol H 2 O m −2 s −1 bar −1 ) has the same meaning as g 0 BB setting a lower bound of the conductance value. g 1 MED is similarly a sensitivity parameter, but in different units ( √ kPa) and scale.
In both models, C s is CO 2 partial pressure at the leaf surface after taking account of boundary layer resistance to CO 2 (r b c ) under a given A n .
h s in Ball-Berry model was obtained by solving an equation connecting diffusion pathways of boundary layer and stomatal interface. h a is relative humidity in the air. Humidity inside the intercellular space was assumed fully saturated. The resultant Equation (A30) is quadratic since g s BB has h s in its component as defined in Equation (A26).
D s in Medlyn model was obtained in a similar way. Water vapor pressure at the leaf surface (w s ), in the air (w a ), and in the intercellular space (w i ) were used instead of relative humidity values. The resultant Equation (A31) is quadratic since g s MED has D s in its component as defined in Equation (A28). In an actual implementation, the equation was solved in terms of √ D s and squared later to cope with a symbolic equation solver internally used by our framework.

. Energy Balance
Leaf temperature (T l ) can be different from air temperature (T a ) due to energy balancing on the leaf surface.
The temperature difference (∆T) is obtained by numerically solving the energy budget equation.
The equation indicates that the net radiation absorbed by the leaf (R n ) should equal to the loss by sensible heat flux (H) and latent heat flux (λE), assuming no heat is stored in the leaf.
The absorbed net radiation (R n ) is composed of shortwave component driven by solar radiation (R sw ) and longwave component via thermal radiation (R lw ). Note that we assumed gas-exchange measured inside a small chamber and the temperature of chamber wall was equivalent to the surrounding air. α s is absorption coefficient of the leaf for solar radiation and is thermal emissivity of the leaf. σ is Stefan-Boltzmann constant.
For the sensible heat flux (H), C p is specific heat of air and g h is the convective heat conductance of the leaf surface.
For the latent heat flux (λE), λ is the latent heat of vaporization for water. E is the transpiration rate calculated with leaf conductance (g v ) and vapor pressure gradient between leaf surface and the air (∆w). e a is ambient vapor pressure of air. e s is saturated vapor pressure at a given temperature (T). The net photosynthetic rate (A n ) depends on mesophyll [CO 2 ] (C m ≈ C i ) through Equations (A2) and (A10). C i depends on stomatal conductance (g s ) through Equation (A6). g s then again depends on A n through Equations (A26) and (A28), forming a cyclic dependency solved by an iterative numerical method. In the meantime, many parameters for A n depends on temperature that the temperature of biochemical reaction site (T l ) should play an important role. Vapor pressure of intercellular space (pi) used for Medlyn stomatal conductance model can also change with T l . The latent heat component of Equation (A33) solved for T l is then driven by leaf conductance (g v ) which depends on g s through Equation (A38).
Therefore, three submodels for photosynthesis, stomatal conductance, and energy balance are interdependent. We used a nested iterative procedure using the bisection method to solve this relation numerically. For initial condition, each variable was given a sensible range of minimum and maximum values. C i was assumed to be at least 0 µbar and lower than two times of C a . ∆T, the difference between leaf temperature (T l ) and air temperature (T a ), was assumed in the range of −10 K and 10 K.
Appendix A.6. Implementation An early version of the model was implemented in C++ and released as a sample of peer-reviewed publication of computer code [58]. Then, the model was translated, calibrated, and visualized with Cropbox framework [50] written in Julia programming language [17]. In this framework, a small component of the model is encapsulated into a structure called system and each system is a collection of variables described in a declarative form which closely resembles mathematical equations defined above and specifications of variables and parameters (Tables A1 and A2). The gas-exchange model was composed of 23 systems including a submodel for stomatal conductance ( Figure A1).
The framework provides a number of unique types of variable declaration that can help hiding complexity of implementation details. For example, h s in Equation (A30) is declared as a solve variable which is automatically expanded and solved in terms of symbolic algebra. ∆T in Equation (A32) in conjunction with Equation (A33) is declared as a bisect variable which would automatically generate necessary code to implement a nested iterative solver, as described in Appendix A.5. An optimal value of ∆T that satisfies R n = H + λE as Equation (A33) would be found out by applying bisection method ( Figure A2). Given Equation (A6) was also implemented as a bisect variable, the framework would analyze dependency between variables and generate proper code for nested structure. Note that all equations described in this section assume implicit units conversion and scaling which is a feature provided by Cropbox framework. For example, in Equation (A23), leaf width W was declared in 'cm' then automatically scaled to match the units of final product in 'm' that an actual code generated for calculating characteristic dimension was d = 0.72 w · 0.01 m cm .  Input variables to the model (Table 1).

Appendix A.7. Workflow
A system implemented on Cropbox framework can be run by simulate() function supplied with an optional configuration. For instance, the coupled gas-exchange model with Medlyn stomatal conductance submodel (GasExchangeMedlyn) was run with a set of default parameters described in Tables 1 and A2 ( Figure A3). simulate() function returns a result of simulation in a tabular data frame which can be further analyzed by various tools including own visualization methods provided by the framework. By default, the result contains values of all numerical variables declared in the controller system.
After testing models with initial parameters, calibration of four parameters, N 0 , s, g 0 MED and g 1 MED , as explained in Section 4.2 can be done by calibrate() function with a training dataset ( Figure A4). In this example, we had two target variables in the dataset, Photo and gs, and wanted to compare them with A_net (A n ) and gs (g s ) estimated by the model. A percentage root mean square error (PRMSE) calculated from the difference between each set of variables was minimized by updating four parameters. The range of each parameter value is specified in parameters option.   Table A2 are included in the baseline configuration c0 used by simulate() function.  (0, 10), N0=(0, 1)), StomataMedlyn  (g0=(0, 1), g1=(0, 10)), ), metric = :prmse, # other options ) Figure A4. Script for calibrating parameters for the gas-exchange model. index option provides a mapping of variable names between a dataset and the model specification. target option indicates two variables, A n and g s , for which fitting errors are minimized. parameters option specifies a range of parameter values to be calibrated. A set of calibrated parameters are stored in a variable named c1.
Plotting model response in multiple ranges of criteria would require aggregation of results from models running multiple times with slightly different configuration. visualize() function takes care of a few common plotting scenarios in a simple interface. All figures in this paper were directly generated from the framework. For example, Figure 2b was generated by running GasExchangeMedlyn model with four levels of RH (RH) from 20 % to 80 % and plotting gs (g s ) against Ci (C i ) in lines ( Figure A5). Each line of response was composed by changing CO2 (C a ) in a range from 10 µbar to 1500 µbar at an interval of 10 µbar.  Figure A5. Script for generating Figure 2b. visualize() function internally runs simulate() and render a plot using the output. Note that a configuration used here is the baseline c0 partially overridden by the calibrated c1. x and y options specify variables used for two axes of the plot. xstep option provides a range of value to compose a single curve, indicating atmospheric CO 2 concentration (C a ) is ranged from 10 µbar to 1500 µbar. group option provides a group of treatments to render multiple curves, indicating four levels of relative humidity (RH) are used. xlim option sets the visible range of C i limited from 0 µbar to 600 µbar and ylim option sets the range of g s from 0 to 1 mol H 2 O m −2 s −1 bar −1 . kind option for plot type can be either :line or :scatter.
A Jupyter notebook containing source code of the model with calibration datasets and scripts for producing figures presented in this paper is available at https://github.com/cropbox/plants2020.

Appendix B. Calibration of SPAD Measurements
Our experimental dataset used in calibration had SPAD measurements as a proxy to leaf nitrogen content. A small number of samples were separately collected from direct measurement of leaf nitrogen content and then used to derive a relationship between SPAD measurement and leaf nitrogen content. A quadratic equation in the form of N = ax 2 + bx + c was fitted for leaf nitrogen content (N) and SPAD measurement (x = SPAD) where a = 0.0004, b = 0.012, and c = 0 with R 2 = 0.92 ( Figure A6). Leaf nitrogen (g m⁻²) Figure A6. Fitting SPAD measurement to leaf nitrogen content (N). N was fitted by a quadratic equation N = 0.0004x 2 + 0.012x for 18 samples of SPAD measurement (x = SPAD) with R 2 = 0.92.

Appendix C. Sensitivity of Nitrogen Parameters
The sensitivity of parameters related to nitrogen dependence was obtained by calculating net photosynthesis rate (A n ) with a range of values for the parameter in question, while keeping other parameters unchanged ( Figure A7). For reference, N 0 was 0.371 and 0.315, and s was 4.470 and 3.912, for BB and MED, respectively. Baseline leaf nitrogen content (N 0 ) barely had an impact on A n with its range from 0.2 g m −2 to 0.5 g m −2 . Sensitivity of nitrogen response curve (s) did not make a large difference either within the range of calibrated parameter values. For example, two sets of N 0 and s mentioned above resulted into less than 1 % difference in A n at 400 µmol m −2 s −1 of ambient CO 2 level. Therefore we assumed the difference between nitrogen parameter values we obtained from calibration for BB and MED were negligible and safe to pool them together for one unified set of parameters. The averaged out values for N 0 and s used in the rest of simulation was 0.343 and 4.191 (Table A2). An (μmol m⁻² s⁻¹) (b) steepness of nitrogen response (s) Figure A7. Sensitivity of parameters related to nitrogen dependence. The coupled gas-exchange model with Medlyn stomatal conductance model was used for testing. Air temperature (T a ) was 32 • C, irradiance (I) was 2000 µmol quanta m −2 s −1 , and relative humidity (RH) was 66 %.

Appendix D. Parameters for Transgenic Plants
A few parameters were derived to replicate a wide range of temperature response of transgenic plants with reduced amounts of Rubisco reported by Kubien et al. [39]. We mostly used the same values for variables and parameters identified in the paper, except J max 25 , which was not mentioned but needed adjustments to replicate original curves (Table A3). Measured values of A n , g s , and C i under three treatments were digitized from their Figure 1. WT indicates wild type, while AR1 and AR2 refers to mutants treated with anti-Rubisco antibody to have Rubisco content reduced to 49 % and 32 %, respectively. Activation energy for V pmax J max 25 120 µmol electrons m −2 s −1 Maximum rate of electron transport at 25 • C V cmax 25 15.54 µmol CO 2 m −2 s −1 Maximum rate of Rubisco carboxylation at 25 • C V pmax 25 146.9 µmol CO 2 m −2 s −1 Maximum rate of PEP carboxylation at 25 • C Minimum g s observed in the dataset. Scaled to the ratio of maximum A n for each treatment to make up missing information.