Kinetic Modeling of Sunflower Grain Filling and Fatty Acid Biosynthesis.

Grain growth and oil biosynthesis are complex processes that involve various enzymes placed in different sub-cellular compartments of the grain. In order to understand the mechanisms controlling grain weight and composition, we need mathematical models capable of simulating the dynamic behavior of the main components of the grain during the grain filling stage. In this paper, we present a non-structured mechanistic kinetic model developed for sunflower grains. The model was first calibrated for sunflower hybrid ACA855. The calibrated model was able to predict the theoretical amount of carbohydrate equivalents allocated to the grain, grain growth and the dynamics of the oil and non-oil fraction, while considering maintenance requirements and leaf senescence. Incorporating into the model the serial-parallel nature of fatty acid biosynthesis permitted a good representation of the kinetics of palmitic, stearic, oleic, and linoleic acids production. A sensitivity analysis showed that the relative influence of input parameters changed along grain development. Grain growth was mostly affected by the specific growth parameter (μ′) while fatty acid composition strongly depended on their own maximum specific rate parameters. The model was successfully applied to two additional hybrids (MG2 and DK3820). The proposed model can be the first building block toward the development of a more sophisticated model, capable of predicting the effects of environmental conditions on grain weight and composition, in a comprehensive and quantitative way.


INTRODUCTION
The weight and composition of oilseed grains at harvest are complex traits that depend on the dynamics of many processes occurring earlier at both the plant and organ levels. Their response to the genotype and the environment results from several linked processes controlled at different levels of organization, from sub-cellular to crop (Martre et al., 2011). The sensitivity of these traits to multiple factors changes during grain development (Aguirrezábal et al., 2003;Rondanini et al., 2003;Echarte et al., 2013). Therefore, understanding how grain weight and composition are determined demands a deep and integrated knowledge of grain filling dynamics.
Most of the photoassimilates supplied to the sunflower grains during their filling period are contemporaneously synthesized by the leaves and thus, leaves are considered the main source of substrate for grain growth (Hall et al., 1990;López Pereira et al., 2008;Echarte et al., 2012). Since grains are the main sink of photoassimilates during this period, changes in assimilate production at the source can be interpreted as changes in carbon availability in the developing grains (Hall et al., 1995), sucrose being the major carbohydrate and the main phloem-transported sugar in sunflower plants (Alkio et al., 2002). Once in the grain, part of this carbon is directed to the pool of acetyl-CoA, the precursor of fatty acids, which are the main components of sunflower oil.
The biosynthesis of fatty acids involves several enzymes placed in different sub-cellular compartments of the grains (Garcés and Mancha, 1991;Gray and Kekwick, 1996;Harwood, 1996). The enzyme acetyl-CoA carboxylase (ACCase) catalyzes the first reaction of this pathway. After successive elongation reactions, the main products of the intraplastidial de novo fatty acid biosynthesis are palmitoyl-ACP (Pp), stearoyl-ACP (Sp), and oleoyl-ACP (Op). These acyl-ACPs are hydrolyzed to free fatty acids, activated to the corresponding acyl-CoAs, and exported to the cytosol to be incorporated into glycerolipids. Finally, oleic acid can be transformed into linoleic acid by the action of oleoylphosphatidylcholine desaturase (FAD2), an enzyme located in the endoplasmic reticulum (E.R), in a step that represents a key point in the regulation of fatty acid composition (Garcés and Mancha, 1991).
Progress in understanding and modeling grain weight and composition dynamics can be achieved with two different main approaches: empirical or mechanistic. While empirical models typically describe the data but do not explain them, mechanistic models are reductionist and explain data based on knowledge of processes at the lower levels of biological organization (Loomis et al., 1979;Thornley and Johnson, 1990). Mechanistic models of physiological processes are often based on biochemical principles such as enzyme kinetics and reaction stoichiometry (Amthor, 2000). For instance, a biochemical model of fatty acids biosynthesis has been proposed by Martínez-Force and Garcés (2002). This model is a structured kinetic one based on every known enzymatic step of the pathway. Although useful to understand the fatty acid biosynthetic pathway, this model does not consider the supply of substrate by the mother plant (which changes along grain filling), nor the synthesis of other grain components that contribute to grain weight and composition. The formulations of this kind of model require extensive knowledge of metabolic pathways, often not available in enough detail, and complex computer programming to carry out the calculations involved (Martínez-Force and Garcés, 2002). Furthermore, the concentrations of intermediate compounds and enzymatic activities are often difficult to measure, hindering validation of state variables, and making them less attractive for practical applications (Durruty et al., 2012).
Mathematical models with a more empirical approach have been developed to predict sunflower development, yield, and yield components (Chapman et al., 1993;Steer et al., 1993;Villalobos et al., 1996;Yeatts, 2004). These crop simulation models are useful tools for evaluating different agronomic management strategies (Villalobos et al., 1996). Based on a certain crop physiology background, they often adequately address the crop growth and development and their interaction with the environment. However, most of these models describe grain filling with insufficient detail, fail to take into account processes occurring inside the grain and rarely consider the accumulation of the main grain components. Some published sunflower models predict oil yield and quality but they are mainly based on empirical relationships between many traits and environmental factors and, moreover, the simulation of grain components accumulation is not dynamic (Pereyra-Irujo and Aguirrezábal, 2007;Casadebaig et al., 2011). To the best of our knowledge, a model describing the dynamics of sunflower grain filling and the accumulation of the main grain components in detail (oil and its fatty acid composition) has not been developed so far.
Non-structured models provide a trade-off between the realism of the biological processes and the relative simplicity required by modeling (Tolla et al., 2007). They are useful tools when access to the data is limited or the complexity of reactions in the pathway hinders the modeling process (Steer et al., 1993;Tolla et al., 2007;Durruty et al., 2012). The aim of the present work was to develop a non-structured mechanistic kinetic model of grain growth and oil and fatty acids biosynthesis, as a tool to describe the dynamics of grain filling in a comprehensive and quantitative way. Such a kinetic model can predict the dynamics of weight and components of the sunflower grain, and would contribute to understand the underlying mechanisms and the response of grain weight and composition to different growing conditions.

General
Our model was built on the basis of biochemical reactions engineering principles. In this first version, the effect of temperature in plant development was taken into account by means of thermal time calculations. Other effects of temperature and effects of variations in daily incident radiation were not considered. Water and nutrients available are considered nonlimiting for all processes. The following main assumptions were made: -the temporal variable t is expressed as thermal time after flowering (accumulated degree days); -grain is considered as a control volume that changes as the grain grows; -grain filling is treated as a fed batch system where the grain grows in batch mode with an external carbon source; -the external carbon source changes with time as a consequence of leaf senescence; -once in the grain, carbon has two possible fates: (i) it turns into substrate for growth (i.e., contributes to grain weight), or (ii) it is used for maintenance.

Grain Growth and Maintenance
Logistic functions have been reportedly useful to predict the sigmoid behavior (S-shape) of sunflower grain weight dynamics (Yeatts, 2004). In the present work, a logistic equation driven by time expressed in degree days (1) was fitted to experimental data to simulate grain growth. Equation (1) was first presented by M' Kendrick and Pai (1911).
In this equation, W is the dry weight of an individual grain, W 0 is the initial amount of W at which the model starts to work and W max is the maximum potential grain weight. The parameter µ ′ represents the specific grain growth rate (i.e., the biomass production rate per unit of biomass). The rate at which the grain grows (r W ) can be written as: The grain is not a closed system as it continuously receives the substrate assimilated by the mother plant. Once in the grain, C substrate rapidly reacts to form the different grain components. Thus, in a C mass balance, the C that enters the grain per unit time is equal to the consumption rate (r C ), while the accumulation rate can be neglected. Thus r C represents the rate of carbon allocation to the grains (or feed rate), here expressed as carbohydrate equivalents per degree days (Vertregt and Penning De Vries, 1987;Echarte et al., 2012). However, not all the substrate reaching the grain is used for growth (defined as increase in grain weight). Part of this C is used to provide the energy needed for diverse processes that do not result in a net increase of dry weight (e.g., turnover of structures, activity of transport and movement, maintenance of concentration gradients and defense systems). The C costs of some of these processes are considered in calculations of carbohydrate equivalents (e.g., maintenance of the tools for biosynthesis, Vertregt and Penning De Vries, 1987). The Pirt's maintenance equation (Pirt, 1975) allows us to separate the C consumed by these processes from the C consumed for growth as follows: where Y G and m are the actual growth yield and maintenance coefficients, respectively. The coefficient Y G represents the biochemical efficiency of transformation of glucose into new plant material (Van Iersel and Seymour, 2000) and m, the C expended in processes that do not result in a net increase in grain dry matter. The coefficient m differs from those previously reported (Ploschuk and Hall, 1997;Van Iersel and Seymour, 2000) in that it does not consider all processes related to maintenance respiration (e.g., cell structure maintenance; Vertregt and Penning De Vries, 1987) since they were not included in carbohydrate equivalents calculations. The first term of Equation (3) represents the substrate used for grain growth (r CG ) and the second one the substrate used for maintenance purposes (r Cm ).Considering total grain weight W as the sum of oil (W O ) and non-oil (W NO ) fractions, each production rate can be obtained by defining Y GO and Y GNO as the actual oil and non-oil fraction yield coefficients, respectively (see Supplementary Material).

Solar Radiation Interception
As the plant life progresses, the capacity of the source to feed the grain decreases because of leaf senescence, and then, r C also decreases. The model takes into account this phenomenon by considering the feed rate as a function of the proportion of the photosynthetically active radiation (p PAR ) intercepted by the crop. The value of p PAR is equal to one at the initial time (beginning of grain filling) and decreases as plant leaves senesce. To solve the model, a mathematical expression of p PAR is necessary. The radiation interception has been previously calculated as a function of the leaf area index, in agreement with the Lambert-Beer law, which in turn depends on thermal time (Gardner et al., 1985;Pereyra-Irujo and Aguirrezábal, 2007). p PAR can be expressed as: where k λ is an empirical parameter that considers the Lambert-Beer law constant and the relationship between the leaf area index and thermal time, and Kp is a proportionality constant between p PAR and radiation interception. In the present work, both parameters were obtained by fitting Equation (4) to experimental data.

Grain Filling
Thus, according to the previous section and in order to take into account the contribution of C assimilated by the mother plant to grain filling, the theoretical accumulation of carbohydrate equivalents, i.e., the cumulative amount of substrate that enters the grain expressed by Equation (3) must be modified to: On the same way, the production rate of W must be expressed as: Finally, according to partitioning of C into the oil and non-oil fractions, the production rate of each component can be depicted as: Equations (5-8) are ordinary differential equations (ODE) and must be simultaneously solved to predict the profiles of C and W during grain filling. The production rate and the substrate used for each component growth and maintenance (r WO, r WNO, r CNO , r CO , and r Cm ), can also be individually predicted by coupling and solving their respective differential equations (see Supplementary Material).

Fatty Acids Biosynthesis
A simplified model of fatty acids biosynthesis is shown in Figure 1. Compounds that react in the plastid are grouped in a global variable C F , which represents the precursor for all the fatty acids produced and stored. This is a simplification of the model presented by Martínez-Force and Garcés (2002). Lumping several intermediates into a global variable C F allows considering the serial-parallel nature of the synthesis pathway without the need for a complex segregated model. Furthermore, the global variable C F implicitly considers the dynamic channeling model presented by these authors. Figure 1 shows the proposed kinetic model and the simplified reactions involved in the synthesis of fatty acids. In previous sections, substrate consumption was divided into maintenance and growth, and the grain growth was fractionated into oil and non-oil components. Since the oil fraction is composed of several fatty acids and intermediates, the rate depicted in Equation (7) pertains to the first step of fatty acids biosynthesis, i.e., the production of intermediary C F . The intermediate C F accumulated is equal to the difference between its production from C allocated to the grain and its consumption to produce P, S, and O, and storage. In agreement with the proposed model, the following net intermediate production rate results: Assuming the production and active transport of P, S, and O follows a specific Michaelis-Menten kinetics (i.e., they depend on grain weight, the heavier the grain, the higher the rate of synthesis of fatty acids), the production rate can be written as: where v max and Ks are the maximum specific rate and half saturation constants, respectively. The subscript i represents P, S, or O. The fatty acid production rate is expressed here as a function of intermediate concentration (i.e., the amount of C F per unit weight). As stated above, linoleic acid is produced from oleic acid inside the endoplasmic reticulum. Then: Finally, the kinetics of C F and the different fatty acids in the oleosome can be calculated via their respective mass balances.
Since W grows with C F , P, S, O, and L, the ODE Block (12-16) and differential Equation (6) must be solved simultaneously in order to predict the biosynthesis of fatty acids during grain filling.

Experimental Data
Sunflower (Helianthus annuus L.) was grown in the field at Balcarce Experimental Station (Unidad Integrada Balcarce INTA-FCA; 37 • S, 58 • W), in the Buenos Aires Province, Argentina. The model was initially calibrated with hybrid ACA885, and later evaluated with experimental data from hybrids MG2 and DK3820. The soil was a Typic Argiudoll. Experiments were performed during growing seasons 2007-2008 and 2012-2013. Each experimental unit consisted of six rows 6 m long spaced at 0.7 m. Plant population density at sowing was 6.5 plants m −2 . The crops were grown under optimal nutrient and water conditions. Soil fertility in all experiments was adequate to attain maximum yields for sunflower crops grown under nonlimiting water conditions-yield >5000 kg ha −1 (Sosa et al., 1999;Andrade et al., 2000). Pests, diseases and weeds were successfully controlled. Flowering of a plant was defined by the appearance of stamens in all florets from the outer whorl of the capitulum-R5.1 stage, (Schneiter and Miller, 1981).

Sampling and Chemical Analysis
Sampling and chemical analysis were performed as in Echarte et al. (2013). Briefly, 12 grains of rows 6-8 were excised from the same plant as long as the total removal did not exceed 5% of the average final capitulum grain number. The number of sampling dates varied between 8 and 12, depending on the experiment. Grains were oven-dried at 60 • C and weighed. Lipids were extracted with 5 ml of hexane:isopropanol (7:2, v/v) and 2.5 ml Na 2 SO 4 (67 g l −1 ) in the presence of 0.2 ml of 1, 2, 3 triheptadecanoyl-glycerol (50 mg ml −1 ) as internal standard. The lipidic phases from the extracted samples were evaporated to dryness under a nitrogen stream. The residue was dissolved in 0.5 ml of hexane and incubated for 1 h at 80 • C in the presence of the methylation mixture methanol:toluene:H 2 SO 4 (88:10:2) and 1 ml of heptane. After samples had cooled down to room temperature, the upper phase containing fatty acids methyl esters was separated. The fatty acid composition of the extracts was determined by gas chromatography (GLC) with a Shimadzu GC-2014 chromatograph (Kyoto, Japan). The fatty acid content and total lipids extracted were calculated with the internal standard method. The oil content was assumed equal to the total extracted lipids given that they reportedly represent more than 96% of the oil (Robertson et al., 1978).

Measurements
Global daily incident radiation was measured with pyranometers (LI-200SB, LI-COR, Lincoln, NE) from a meteorological station located ∼400 m away from the experimental units. The proportion of photosynthetically active radiation (pPAR) intercepted by the crop at noon (±1 h) was calculated according to Gallo and Daughtry (1986) as (1 − Rb/Ro), where Rb is the radiation measured below the oldest green leaf, and Ro is the radiation measured above the canopy. Rb was measured weekly with a line quantum sensor (LI-191SB, LI-COR, Lincoln, NE, USA) positioned across the rows (the length of the sensor was modified according to the distance between rows, 0.7 m). Three measurements were taken per plot. Air temperature was measured using shielded thermistors (Cavadevices, Buenos Aires, Argentina) next to the capitulum every 60 s and averaged hourly. Measurements began after flowering and finished at physiological maturity and were recorded by data loggers (Cavadevices, Buenos Aires, Argentina). The amount of assimilates effectively allocated to the grains (C) was assumed to be represented by carbon equivalents for grain biomass production (Vertregt and Penning De Vries, 1987). For this, carbon and nitrogen in the grains were determined with a TruSpec CN equipment (Leco Corporation, St. Joseph, MI), and the ash content was measured according to AOAC recommendations (Aoac, 1990). Carbohydrate equivalents for grain biomass production were calculated as described by Vertregt and Penning De Vries (1987).

Thermal Time
The temporal variable t is expressed as cumulative degree days, with the aim of expressing time and rates in a temperature compensated way to make temporal effects independent of temperature fluctuations (Kiniry et al., 1992;Parent and Tardieu, 2012). Cumulative degree days are calculated from daily data for mean temperature (Tm) and a base temperature (Tb) of 6 • C (Kiniry et al., 1992) as follows: Values of C F were calculated for every "j" thermal time from experimental data using a box mass balance Equation (18). Y WO/C represents the apparent oil yield coefficient, it was obtained from the linear regression of oil vs. C experimental data.

Actual Yield Coefficients
The actual growth yield coefficient (Y G ) and the maintenance coefficient (m) were calculated by linear fitting of Pirt's Equation (Pirt, 1975) to experimental data: where Y is the instantaneous growth yield (Pirt, 1975). Values of Y and µ were calculated for every "j" time interval as central differences of the experimental data in agreement with Equations (20) and (21), respectively.
Subsequently, the actual non-oil fraction yield coefficient (Y GNO ) and actual oil fraction yield coefficient (Y GO ) were determined using a general non-linear regression.

Kinetic Parameters
Growth kinetic parameters W 0 , W Max , and µ ′ were calculated by non-linear fitting of Equation (1). The kinetic parameters for fatty acids production v max and Ks were obtained by non-linear fitting of Equations (10) and (11). The values of production rates were obtained by finite differences (central differences) according to Equation (22): where i represent P, S, O, or L. Once all the parameters were obtained, a further step of model optimization was performed using a general non-linear regression framework.

Sensitivity Analysis
A sensitivity analysis was performed as in Villalobos et al. (1996). The effects of ±25% variation in every kinetic parameter of the model on the main variables output (C, W, P, S, O, and L) were analyzed. The sensitive coefficient (SC) was calculated as in Equation (23), being V the output variable and P the kinetic parameter.

Informatics Tools and Statistics
Linear and non-linear regressions were performed with Origin 8.0 R (OriginPro, v. 8.0724; OriginLabCorporation, Northampton, MA 01060, USA). Once the kinetic parameters were obtained, profiles were modeled using a fourth-order Runge-Kutta algorithm coupled to the regression in order to integrate the differential equations simultaneously (MathCad 14.0.0.163, Parametric Technology Corporation). The goodnessof-fit of the model was evaluated using the regression coefficient (R 2 ) and Reduced Chi-square (χ 2 ) test via the Prob (χ 2 > F) with α < 0.05. Significant differences among parameters were evaluated by Student t-test (95% confidence interval).

Grain Growth and Filling
Grains grow with thermal time following a sigmoid curve ( Figure 3A). After fitting Equation (1) to the experimental data, the values obtained for W 0 , W max and µ ′ were 0.4299 ± 0.1634 mg, 34.5396 ± 0.9131 mg, and 0.0145 ± 0.0013 • Cdaf −1 , respectively [R 2 = 0.94547, Prob (χ 2 > F) = 0]. Figure 3B shows the experimental values of C vs. thermal time. The values obtained for Y G and m by fitting Equation (19) to the experimental data were 0.671 ± 0.036 mg mg −1 C and 4.301 × 10 −3 ± 4.28 × 10 −4 mg C mg −1• Cdaf −1 , respectively. The fitted values for Y GNO and Y GO were 0.368 ± 0.021 g g −1 C and 0.303 ± 0.016 g g −1 C , respectively [R 2 = 0.9552, Prob(χ 2 > F) = 0]. Grain weight and cumulative carbohydrate equivalents were simulated by solving simultaneously ODE Block (5-8) with the previously fitted parameters. The values predicted by the model are presented in Figures 3A,B as solid lines. The figures indicate the model successfully describes the experimental grain growth kinetics and it is able to predict the theoretical cumulative amount of carbohydrates that is allocated to the grain during the filling period. The amount of substrate spent in growth (CG), maintenance (Cm), oil (CO), and non-oil (CNO) grain fractions were also predicted and plotted as a function of thermal time (dashed lines in Figure 3B). Different dynamics of carbon investment were observed: at first, most of the substrate is used to produce grain mass, later, the maintenance requirements increase and thus less substrate is used for growth.
Values of W NO and W O were also obtained by solving ODE Block (5-8) (Figures 3C,D). The simulations seem to overestimate the experimental W O -values before 200 • Cdaf.
Grain weight was plotted as a function of carbohydrate equivalents in Figure 4. The solid line shows the values predicted by the model; the instantaneous growth yield is represented by the slope of this line. It can be observed that the model acceptably simulates the non-linear behavior of experimental data. At the beginning of grain filling, when grains are small, they present lower maintenance requirements, which results in a higher apparent growth yield. As grains grow in size, more substrate is used for maintenance and the apparent yield coefficient decreases. Table 1 shows the kinetic parameters obtained from fitting Equation (10) and (11). The values predicted by the model are presented in Figure 5 as solid lines. The model successfully describes the experimental behavior of all fatty acids. The rates of production of every fatty acid increase early during grain filling together with grain weight, later they decrease until they completely stop. The model predicts the amount of oleic acid increases at the early stages of grain filling, reaches a maximum at 400 • Cdaf and then decreases. Linoleic acid production follows an end-product saturating specific kinetics, like palmitic, and stearic acid, but smoothed and delayed by a reversible reaction and by the fact that linoleic acid is the final product of three serial steps (C→C F →O→L).

Fatty Acids Biosynthesis
Once the profile of each fatty acid was obtained, the concentration of oil was calculated as the sum of all fatty acids and plotted in Figure 5E as a solid line. The dotted line in the figure represents the oil grain weight predicted in a first approach (Figure 3D), where a single step production process was considered (see Section Grain Growth and Filling). A better performance of the model was achieved when the serial nature of fatty acids biosynthesis was considered, especially at short times (<300 • Cdaf) when C F is accumulated.

Sensitivity Analysis
A sensitivity analysis was performed on the proposed model. The effect of ±25% variation in every kinetic input parameter on the model output was observed at three different times during grain filling. The results, presented in Table 2, show that the relative influence of the parameters on both grain weight and composition varies along the grain development. Specific growth (µ ′ ) has significant influence on all the variables at the beginning of grain filling (SC > 0.5) but its effect decays later on (SC < 0.5), when C is used in maintenance processes. On the other hand, W 0 , W max , Y G , and m exert low influence (SC < 0.3) on most of the variables. However, palmitic, stearic, and linoleic acids were sensitive to grain growth and oil yield (Y G and Y GO ), showing the influence of the partitioning of carbon to the oil or non-oil fraction on fatty acid composition. Saturated fatty acids were sensitive to their own maximum specific rate of synthesis (v maxS and v maxP ) and to the oleic acid maximum specific rate (v maxO ), while oleic acid was mostly influenced by parameters driving the synthesis of linoleic acid.

Model Extrapolation to Different Hybrids
The ability of the model to predict the behavior of two independent hybrids that were not used for model calibration (MG2 and DK3820) was evaluated. Parameter values and their comparison among hybrids are provided as Supplementary Material (Table SM1). In order to define the smallest set of values necessary to obtain an appropriate simulation, the model was run by: (i) freely fitting all the parameters, or (ii) fitting the minimum amount of parameters that would ensure a good predictive quality [Prob(χ2 > F) < 0.01]. For this, parameters of low sensitivity (W 0 , W max , Y G , m, Ks i ) and µ ′ (sensitive only at the beginning of grain filling and low genetic variability) were fixed and v maxi parameters were refitted. In (ii), the ACA885 parameter values were used instead of each hybrid's own parameters (bold values in Table 3). Table 3 summarizes the kinetic parameters obtained by fitting the model to the experimental values measured for all hybrids.
Growth parameter values are similar for all hybrids, with the exception of Y G and Y GO of MG2. Parameter µ ′ showed high sensitivity early during grain filling but did not differ among hybrids, while v maxO showed both low sensitivity and low genetic variability. When the set of parameters (i) was used, the model successfully predicted the kinetics of every trait explored in this research (data not shown). Grain weight, theoretical accumulated carbohydrates, and fatty acid composition dynamics simulated by fitting the minimal amount of parameters (case ii in Table 3) are shown in Figure 6. In this case, simulated values of C and W adequately described the experimental behavior of both MG2 and DK3820 hybrids (Figures 6A,B). On the other hand, the extrapolation of all fatty acids biosynthesis parameters was not possible due to the high sensitivity of P, S, O, and L to v maxi . The model adequately predicted fatty acids dynamics during grain filling when five out of the 16 parameters were refitted.

DISCUSSION
In the present work, a non-structured mechanistic kinetic model of grain growth and oil and fatty acids biosynthesis has been developed. By setting initial conditions -W 0 for grain weightand calculating carbon assimilated by leaves and allocated to the grains as the substrate, the oil, and non-oil weight and oil composition dynamics have been successfully simulated for different sunflower hybrids. To the best of our knowledge, a model with the ability of describing the grain filling dynamics in such detail has not been previously developed for sunflower, or any other crop species.

Carbon Partitioning to Growth and Maintenance
The model considered that the carbon substrate was destined to both growth and maintenance. The amount of substrate that was actually transformed into grain biomass (Y G ) was in the range of previously reported values for conversion efficiency (0.6 to 0.8 mg mg −1 Mccree, 1982;Van Iersel and Seymour, 2000). The maintenance coefficient (m) was as well in the range of reported values [0.003 to 0.050 mg mg −1 d −1 (Hesketh et al., 1980)], despite in the present work it did not consider the carbon costs of cell structure maintenance (Penning De Vries et al., 1974). High variability of m has been associated to the dependence of this coefficient on the age of the plant and the environmental conditions during grain filling (Van Iersel and Seymour, 2000).
The results show that if maintenance processes were negligible, 45% of the carbohydrates destined for growth would be transformed into oil. When maintenance processes were considered, 45% of the carbohydrates allocated were destined to grain growth and only 40% of them were converted into oil. Therefore, the results of the model indicate that maintenance processes not only reduce the grain growth, but also the selectivity to oil. In light of these findings, further research might help to understand the physiological processes underlying the relationship between maintenance and grain composition.
Carbon partitioning to maintenance or growth changed with ontogeny. As the grain grew, the substrate destined to maintenance increased ( Figure 3C) and the apparent yield coefficient decreased, in agreement with Pirt's law, (Equation 3;Pirt, 1975). Van Iersel and Seymour (2000) found that r Cm depends on both the age of the plant and the biomass dry weight. A similar behavior was found when analyzing the consumption rate of substrate for maintenance (r Cm ) or grain growth (r CG ) as a function of thermal time (Figure 7A) or grain weight (Figure 7B; see Supplementary Material for r Cm    Positive and negative change refers to the fitted kinetic parameter ±25% variation, respectively. Frontiers in Plant Science | www.frontiersin.org and r CG calculation). The value of r CG presents a maximum at 300 • Cdaf, in concordance with the inflection point of the sigmoid growth curve ( Figure 3A). The value of r Cm reaches its maximum later and at higher grain weight than r CG (Figure 7B). The earlier decrease of r CG indicates the system is more selective to maintenance when the grain is bigger. In this sense, Van Iersel and Seymour (2000) propose that younger plants do not show substrate limitations and maintenance increases as the grain grows. The amount of substrate consumed for growth (r CG ) increases with W due to higher carbon use efficiency as the plants become bigger. As the plant life progresses (later in the plant cycle and higher W), the substrate available diminishes together with the substrate destined for both, growth and maintenance. According with Figure 1A when r C falls due to p par effect, both r Cm and r CG fall because they are also limited by carbon allocation (Equation 5).

Simulation of Grain Weight and Composition Dynamics
Grain growth (total weight, oil, and non-oil components) followed sigmoid functions with time, in agreement with many reports in the literature (Aguirrezábal et al., 2003;Mantese et al., 2006;Rondanini et al., 2007;Echarte et al., 2013). In a previous work, Echarte et al. (2012) reported that the grain weight and oil content linearly increased with the amount of carbon allocated to the grains. In this research, a model with a more mechanistic approach was able to predict the theoretical cumulative amount of carbohydrates that was allocated to the grain during the filling period, and successfully described the grain growth kinetics.
The kinetic parameters of fatty acid biosynthesis were first obtained for sunflower hybrid ACA885. Given a parallel reaction scheme, the system is considered more selective toward the reaction with higher rate. A higher maximum rate of production for oleic acid (v maxO ) than for palmitic and stearic acids made the system more selective toward oleic acid. Similar values of Ks for P, S, and O, which represent the affinity of enzymes involved in their active transport out of the plastid, suggest that these enzymes have similar affinity for the three fatty acids. According to the model predictions, palmitic, and stearic acids followed typical end-product saturating specific kinetics (Figures 5A,B), while oleic acid increased at early stages of grain filling up to a maximum. These predictions are in agreement with previous studies (Martínez-Force et al., 1998;Santonoceto et al., 2003;Echarte et al., 2013). One possible explanation for the behavior of oleic acid is that as grain filling progresses, oleate desaturase activity increases (Gray and Kekwick, 1996), but carbon accumulates in the grain faster than the increment of this catalytic activity. Therefore, between 150 and 300 • Cdaf, oleic acid begins to accumulate. Between 300 and 350 • Cdaf, a high desaturation activity produces a decrease of oleic acid with a concomitant increase of linoleic acid, being more evident when carbon is scarce (Echarte et al., 2013). The accumulation of linoleic acid responds to higher values of v maxL than Vmax−L . Furthermore, v maxL was the highest v max value of all, explaining the high productivity of this compound in agreement with previously reported data (Martínez-Force et al., 1998;Santonoceto et al., 2003;Echarte et al., 2013).
In a first approach, the oil fraction was predicted as if it were produced in a single reaction step (Figure 3D), although it is the final product of a more complex reaction pathway. The low oil fraction values estimated this way suggest an accumulation of intermediate compounds, that were assigned to the non-oil fraction in this first approach. However, when total fatty acids (oil) were estimated as the sum of every fatty acid predicted ( Figure 5E) a better performance was achieved. Thus, the model depicted in Figure 1, which considers the serial nature of fatty acids biosynthesis, was able to better represent the experimental FIGURE 6 | Grain filling dynamics of different hybrids. Simulations were done using parameters set (ii) depicted in Table 3 behavior, especially at short times (<300 • Cdaf) when C F is accumulated.

Finding Key Genotype Parameters of the Model
Complex models are not suitable for the characterization of the dynamics of multiple genotypes, since once the model has been built, finding the kinetic parameters of a new hybrid might be laborious, expensive, and time consuming. Thus, many models rely on a limited number of genetic parameters that appropriately describe one particular genotype behavior, while assuming that the rest of the parameters do not significantly influence the model output for any genotype (Quilot et al., 2005;Makowski et al., 2006). Given a small set of parameters, re-parameterizing many growing models may not require new dynamic measurements, and parameters can be estimated from the final values by optimization methods. The model developed here for sunflower hybrid ACA855 used 16 input parameters to simulate grain weight and component dynamics. Combining the results of sensitivity analysis and the genetic variability of parameters of two other hybrids (MG2 and DK3820), the number of model input parameters was reduced to five. Whether this is sufficient to simulate the behavior of the universe of commercial hybrids should be further tested working with a bigger pool of hybrids than the one explored in this research.
A sensitivity analysis showed that the influence of parameters changed with ontogeny. Although the specific growth rate (µ ′ ) was the most influential grain growth parameter on model output at the beginning of grain filling, similar values among hybrids indicate low genetic variability and thus, a unique value of µ ′ could satisfactorily simulate grain filling dynamics for any hybrid. Although they are parameters of low sensitivity, higher values of conversion efficiencies (Y G and Y GO ) for MG2 suggest this hybrid devotes more substrate to grain weight than the other hybrids. Fitting of parameters related to fatty acids biosynthesis, more specifically the maximum specific rates (v maxP , v maxS , v maxO, v maxL , v max−L ), is needed to re-calibrate the model for every hybrid to adequately predict the fatty acid composition. Traditional commercial hybrids have been improved to obtain maximal grain weight and oil content, but not targeting their fatty acid composition. The latter could have been modified or unintentionally selected when domesticating or breeding other characters (Chapman and Burke, 2012). Hybrids with different potential fatty acid composition have been obtained by mutagenesis (e.g., high oleic hybrids) and oils with certain fatty acid composition receive a prime over the regular price in the market. This model could help to understand the dynamics of oil and fatty acid biosynthesis in sunflower hybrids with modified potential fatty acid composition (high oleic, high stearic, high oleic-high stearic, etc.).

Potential Uses of the Model
The model presented in this paper has been mainly targeted at simulating genetic effects on grain filling and composition dynamics. The identification of key genotypic parameters could guide future research on physiological processes and guide breeding programs. In addition, it is a promising tool to model the effects of biotic and abiotic factors on these dynamics. In this first version, the assimilates availability was estimated based on available data of canopy interception. However, the model could be linked to crop models capable of simulating it by considering different environmental input variables and calculating intermediate phenological or plant structure ones (e.g., Pereyra-Irujo and Aguirrezábal, 2007). In addition, using the present model as a platform and making the necessary modifications, it will be possible to explore the dynamics of other oilseed species (like soybean or rape), where enzymes and pathways are known to significantly differ from those in sunflower.

CONCLUSION
In this work, a kinetic model of sunflower grain filling and fatty acids biosynthesis has been developed. The ability of the model to predict the experimental values was successfully evaluated and validated in different hybrids. The combination of sensitivity analysis and the genetic variability of parameters allowed minimizing the number of input parameters required to appropriately simulate the dynamics of grain filling and component accumulation in different hybrids. The growth model considered a simple effect of carbon source dynamics, maintenance requirements, and a simplified serial-parallel reaction system to describe the fatty acids biosynthetic pathway. The model developed represents a useful tool for future research to evaluate the effects of different factors on grain weight and composition, in a comprehensive and a quantitative way.

AUTHOR CONTRIBUTIONS
In this work the skills from two different knowledge areas have been integrated. The experience of ID on kinetic modeling based on chemical reaction engineering basis was used to modelate experimental data from agricultural field. The experimental data have been obtained in INTA-BAlcarce by ME and LA during several years. ME provided the theoretical base and ID developed and programmed the kinetic model, aiming to represent the metabolic behavior of sunflower grain filling and fatty acids biosynthesis. All the authors wrote the manuscript.