Defining the scope for altering rice leaf anatomy to improve photosynthesis: a modelling approach

Summary Leaf structure plays an important role in photosynthesis. However, the causal relationship and the quantitative importance of any single structural parameter to the overall photosynthetic performance of a leaf remains open to debate. In this paper, we report on a mechanistic model, eLeaf, which successfully captures rice leaf photosynthetic performance under varying environmental conditions of light and CO2. We developed a 3D reaction‐diffusion model for leaf photosynthesis parameterised using a range of imaging data and biochemical measurements from plants grown under ambient and elevated CO2 and then interrogated the model to quantify the importance of these elements. The model successfully captured leaf‐level photosynthetic performance in rice. Photosynthetic metabolism underpinned the majority of the increased carbon assimilation rate observed under elevated CO2 levels, with a range of structural elements making positive and negative contributions. Mesophyll porosity could be varied without any major outcome on photosynthetic performance, providing a theoretical underpinning for experimental data. eLeaf allows quantitative analysis of the influence of morphological and biochemical properties on leaf photosynthesis. The analysis highlights a degree of leaf structural plasticity with respect to photosynthesis of significance in the context of attempts to improve crop photosynthesis.


SI Methods: THE eLEAF MODEL THEORY GUIDE FOR THE eLEAF MODEL
Workflow of eLeaf, and inputs & outputs between modules.

-Module of 3D reconstruction
In order to reconstruct a 3D leaf geometry with exactly the same anatomical features as the measured values (see Table 1), a delicate design of cell shape and cellular structure was developed, as shown in the accompanying diagram of a mesophyll cell with six lobes, where the structure is controlled by parameters a-l and number of lobes of the mesophyll cell.
-'a' controls the thickness of cytosol between chloroplast layer and cell wall 4. Since the measured number of lobe equals to four, the maximum number of lobe allowed for IR64_aCO2 model and IR64_eCO2 model was set to be six. 5. Next, a frame of leaf unit was built based on leaf/mesophyll thickness at minor vein and at bulliform cells, and also shape and size of bundle sheath cells. Then, we started to fill in the region of mesophyll with MCs. (see animation clip -'Match mesophyll porosity').
6. MSCs with different number of lobe were filled into the space in a compact way to reach a leaf porosity value closest to the measured one. Then thickness of the modeled leaf unit was adjusted to match the porosity exactly, considering the laminar structure of rice leaf in the longitudinal direction With the reconstructed 3D leaf, triangle surface meshes were generated for all cells and cellular structures including chloroplast layer and vacuole. Mitochondria were neglected in current ray tracing to reduce the computational complexity. The same ray tracing algorithm as Xiao et al., 2016 was applied here to simulate the light propagation inside leaf and predict light absorptance of each chloroplast in both mesophyll cells and bundle sheath cells.
To compare with the measured physiological data from Licor 6800 (LI-COR Inc., Lincoln, NE, USA), the light source was set to be collimated light with 90% red light (625nm) and 10% blue light (475nm). Special boundary conditions were assigned to the front, back, left, and right surfaces of the leaf unit, so that rays hit these boundaries would be reflected back into the leaf unit. In another word, only rays hit the upper or lower boundary would be counted as total reflectance and transmittance respectively.
Reflection and refraction of light on the interface of two substances were calculated based on Fresnel equations. Meanwhile, light absorption inside certain substance was calculated from Beers-Lambert's law. Same refractive indexes and absorption coefficients of cell walls, cytosol, chloroplasts, and vacuoles as Xiao et al., 2016 were adopted. Absorption coefficient of chloroplast is also linear to the chlorophyll concentration. Here by assuming the chlorophyll concentration distributes uniformly across the leaf, we can calculate it from the measured chlorophyll content per leaf area and the information about chloroplast volume based on the reconstructed 3D leaf geometry.

-Module of CO2 reaction-diffusion and photosynthetic metabolism
A constant CO2 concentration ([CO2]) was set in the substomatal cavity in the 3D geometry as part of the boundary condition. For the purposes of eLeaf, the rest of the epidermis was assumed to be impermeable to CO2. In the intercellular air space, CO2 molecules move following the diffusion equation. On the outer surface of the cell wall, [CO2] in the gaseous phase was converted to [CO2] in the liquid phase following Henry's Law.
The internal surface representing the cell wall and plasma membrane was modelled as a thin diffusion barrier with a given permeability for CO2, (S1) where n is the normal vector of the surface, therefore the left-hand side of the equation is the diffusive flux. PCO2 (m s -1 ) is the permeability of the boundary to CO2, and C1 and C2 (mol m -3 ) are the concentration on either side of the surface. Internal boundaries representing the chloroplast envelope, mitochondria envelope and tonoplast were modelled in the same way except that the cell wall and membrane were set to be impermeable to HCO3 -.
In the compartments of the cytosol, chloroplast, mitochondria and vacuole, reactiondiffusion processes of CO2 were described by a general equation: where Dc (m 2 s -1 ) is the liquid-phase diffusion coefficient of CO2 in water, rf,i is a dimensionless factor representing the decrease of the diffusion coefficient relative to free diffusion in water in different compartments.
is the Laplace operator which equals . While on the right-hand side of the equation, f is volumetric carboxylation rate (mol m -3 s -1 ), h is hydration rate from CO2 to HCO3catalyzed by CA, and rd is volumetric respiration rate, and rp is volumetric photo-respiration rate. In addition, these terms are distributed differently in each compartment, for example, in the cytosol f = rd = rp = 0, in the chloroplast rd = rp = 0, and in mitochondria f = 0. The reaction-diffusion processes of HCO3-in the cytosol, chloroplast, mitochondria and vacuole were described by: where Db (m 2 s -1 ) is the liquid phase diffusion coefficient for bicarbonate, and B (mol m -3 ) is the bicarbonate concentration.
The volumetric carboxylation rate f is calculated based on the FvCB model: Where fc (mol m -3 s -1 ) is volumetric Rubisco limited carboxylation rate, and fj is (mol m -3 s -1 ) volumetric RuBP-regeneration limited carboxylation rate. kc (s -1 ) is Rubisco turnover rate, and Xc (mol m -3 ) is the Rubisco concentration. C (mol m -3 ) is the CO2 concentration, Km The volumetric photorespiration rate rp in the mitochondria was calculated as integration over the chloroplast volume and assigned uniformly to mitochondria by dividing by the mitochondria volume (Vm) in that mesophyll cell (Eqn S5). (S5) The volumetric respiration rate rd, was calculated from the respiration rate per leaf area (Rd) by dividing by the mitochondria volume (Eqn S6). (S6) The hydration rate h (mol m -3 s -1 ) was approximated by Eqn S7.
where ka (s -1 ) is the turnover rate of CA, Xa (mol m -3 ) is the concentration of CA, H (mol m -3 ) is the proton concentration, Ka (mol m -3 ) and KHCO3 (mol m -3 ) are the Michaelis-Menten constants of hydration and dehydration, and Keq (mol m -3 ) is the equilibrium constant.

-Output An and ΦPSII based on the solution from solver
A finite element method was applied to solve the reaction-diffusion system, from which the steady-state CO2 concentration in 3D space was predicted. By integration and summation of the net photosynthesis rate in each cell, it's easy to obtain An for the whole leaf.
The calculation of ΦPSII in eLeaf is modified from the method in Evans 2009, which is used to calculate ΦPSII in a spinach model consisting of 17 paradermal layers.
The light absorbed by PSII in each cell driving linear electron transfer chain, where I is the incident irradiance, ab (i) is the light absorptance of i th cell, β is the proportion of absorbed light partitioned to PSII.
Then , i.e. the light limited rate of PSII electron transfer rate of i th cell, is calculated by . Y(II)LL is the conversion efficient of photosystem II from absorbed photons into e -.
For cells whose photosynthesis is RuBP regeneration-limited, the electron transport rate is calculated as If the photosynthesis of i th cell is Rubisco limited, then its is calculated from (S10) The photochemical efficiency of i th cell, , is therefore (S11) The photochemical efficiency of the leaf, , is measured from the fluorescence under actinic light, Fs, and the fluorescence under a saturating pulse, (Genty et al., 1989). (S12) In a similar way to Evans 2009, Fs equals to the sum of Fs (i) for each cell, while where c is an assumed constant.
equals the sum of for each cell, while is approximated by Eqn. S13, Eqn. S14 and Eqn. S12 combine into (S15) and the assumed c is actually canceled out in the final expression.
In this way, with predicted ab (i) and from eLeaf together with s, Y(II)LL, , and , ΦPSII under different conditions can be calculated.

GENETIC ALGORITHM FOR PARAMETER ESTIMATION
A genetic algorithm (GA) is a metaheuristic search algorithm inspired by the process of natural selection. Five phases are included in a genetic algorithm, 1) initialization of population; 2) evaluate fitness; 3) selection; 4) crossover; 5) mutation, as summarized in the Workflow below.
Workflow of a genetic algorithm. [CA] and wall permeability).
[CA] and wall permeability are assumed to be equal between the IR64-aCO2 and IR64-eCO2 models, so together we have 12 parameters to fit for each individual.
A generation of GA is one iteration of the algorithm from 'evaluate fitness' to generate a new population via 'selection', 'crossover' and 'mutation'. Here we set 12*6=72 as the population size of each generation. The recommended population size is 4-fold to 10-fold larger than the number of variables being fitted.
After setting the upper and lower boundary of each variable, an 'initial population' was generated randomly in the parameter space. Here, to increase the efficiency of searching at the beginning, we also suggest adding an 'individual' obtained by fitting the physiological data under the classic FvCB model.

-Evaluate fitness
Fitness of each 'individual' is the reciprocal of the least squared residual of simulated AN & ΦPSII with experimental data. Considering AN varies between 0 to 50, while ΦPSII is always under 0.9. Direct calculation of the least squared residual will over-estimate the weight of AN. Therefore each response curve (light and CO2 under low or normal oxygen) was normalized to its maximum value before the calculation of the least squared residual.
The closer the simulation is to the experimental measurement, the higher the fitness is. -Selection Ten individuals are selected as parents to build the new generation of population. Normally, a genetic algorithm just take a random strategy during the selection. The probability of choosing an individual as one of the parents is proportional to its fitness. Here we made the probability of choosing one of the top 10 individuals with the largest fitness 3-fold higher, which seems accelerate the convergence of our parameter estimation. 3) [(a1+b1)/2,…, (a12+b12)/2] These three possible new 'individuals' are also randomly selected during the 'crossover'.

-Mutation
After the 'crossover' step, we have a new population of 12*6 individuals. Then a step of 'mutation' is introduced to ensure that it is theoretically possible for the GA to search the whole parameter space. A small mutation rate (5%) of each variable in each 'individual' is implemented. The selected "loci" of mutation will be assigned with a new value randomly from its range. • To repeat the eLeaf IR64-eCO2 model, type § run_e_leaf_v1_2_5(4, [1 1 1 1 1 1 1 1 1]) • The first input number "4" here configures an operating mode, meaning eLeaf will run the simulation under the same light and Ci conditions as the experimental data.
The second input is a vector configuring values assigned to each category of parameters (8 categories in all). 0 means the model will use values from the IR64-AC model, 1 means the model will use values from the IR64-EC model.

Inputs and Outputs of eLeaf
Inputs in Table 1 and Table 2 are assigned in e_geo_parainput_v1_2_5_a4tfit.m Experiment data are in 2.5.eleaf_fvcb_fit/export_meas4comsol_new_selected.m.
We measured the CO2 response and light response under ambient oxygen, correspondingly there are two studies for eLeaf simulation: • Study 1 is the light response curve.
• Study 2 is the CO2 response curve.
Outputs will be in folder eleaf/2.5.eleaf_fvcb_fit/, including • a mph file (COMSOL model file) named "eleaf_fvcb_prswp_CKIR64_a_phipsii.mph" or "eleaf_fvcb_prswp_HCIR64_a_phipsii.mph", which can be opened in COMSOL GUI. This is a file of the model with steady state solutions of both studies.
• "a_std1_prswp_CKIR64.txt" and "a_std2_prswp_CKIR64.txt" are exported files from solutions in the mph file. Predicted leaf photosynthesis rate under different conditions are extracted and recorded.
• "phipsii_std1_prswp_CKIR64.txt" and "phipsii_std2_prswp_CKIR64.txt" are exported files from solutions in the mph file. Predicted carboxylation rate and Cc of each cell under different conditions are extracted and recorded. ΦPSII of the leaf can be calculated based on these files.

Automatic 3D reconstruction module
The current version of eLeaf (v1.2.5) has only been tested with two sets of parameters from IR64-aCO2 and IR64-eCO2 plants. Therefore the author can only ensure that within these ranges, the eLeaf model can run successfully. Extreme input values such as a very small intercellular air space, or a very big or small plastid volume may lead to warnings in the MATLAB command window or even crash the model.

Light propagation module
Information of the boundaries for the ray tracing simulation needs to be passed from the module of 3D reconstruction. For each new 3D geometry, a new 2.e_raytracing/Defs.h file will be generated correspondingly. The triangle meshed surfaces of the reconstructed 3D geometry are represented in the ply format, which are generated by the function  For all data categories, n =5 except for category F3 and F7, n = 15. Relative difference was calculated as (eCO2-aCO2)/aCO2 (%)). BS= bundle sheath; MC = mesophyll cell. Smes = exposed mesophyll cell area (%). The eLeaf model consists of a series of virtual cells bounded by parameters of leaf and mesophyll thickness and interveinal distance set by data obtained from hand-sections of rice. Bundle sheath size and plastid volume, as well as vein size obtained from TEM images are used to set these parameters, with mesophyll cell size and lobing data obtained from confocal light microscopy. Mesophyll cell wall thickness, plastid volume and exposed mesophyll cell area are obtained from TEM images, with mesophyll cell shape and separation adjusted to match overall mesophyll porosity calculated by microCT imaging of rice leaves