Simulating organ biomass variability and carbohydrate distribution in perennial fruit crops: a comparison between the common assimilate pool and phloem carbohydrate transport models

Variability in fruit quality greatly impedes the profitability of an orchard. Modelling can help find the causes of quality variability. However, studies suggest that the common assimilate pool model is inadequate in terms of describing variability in organ biomass. The aim of the current study was to compare the performances of the common assimilate pool and phloem carbohydrate transport models in simulating phloem carbohydrate concentration and organ biomass variability within the whole-plant functional–structural grapevine ( Vitis vinifera ) model that we developed previously. A statistical approach was developed for calibrating the model with a detailed potted experiment that entails three levels of leaf area per vine during the fruit ripening period. Global sensitivity analysis illustrated that carbohydrate allocation changed with the amount of leaf area as well as the limiting factors for organ biomass development. Under a homogeneous canopy architecture where all grape bunches were equally close to the carbohydrate sources, the common assimilate pool and phloem transport models produced very similar results. However, under a heterogeneous canopy architecture with variable distance between bunches and carbohydrate sources, the coefficient of variation for fruit biomass rose from 0.01 to 0.17 as crop load increased. These results indicate that carbohydrate allocation to fruits is affected by both the size of crop load and fruit distribution, which is not adequately described by the common assimilate pool model. The new grapevine model can also simulate dynamic canopy growth and be adapted to help optimize canopy architecture and quality variability of other perennial fruit crops. Results were obtained from the coupled phloem/xylem transport (CT), and no significant difference was found on the percentage of carbon loading and unloading by the grouped organ types between the common assimilate pool (CP) model and CT model in this


IN TROD UCTION
The majority of non-structural carbohydrate (NSC) used by vascular plants is not used where it is fixed or stored, i.e. the source organ, but is transported to other metabolically active organs, i.e. sink. The transport and distribution of carbohydrates are essential for plant survival, vegetative growth and reproductive development (Savage et al. 2016). The allocation of carbohydrates between the vegetative and reproductive organs in perennial fruiting crops is crucial for reaching and maintaining high productivity in terms of fruit yield. Furthermore, carbohydrate distribution to the fruit contributes to product quality. In terms of viticulture, vineyard management inputs, for example, leaf removal, fruit thinning and winter pruning, are commonly conducted with the goal of manipulating the carbohydrate allocation between the different plant organs, through attempting to optimize the ratio between the net source and sink organ sizes during the growing season (Holzapfel et al. 2010;Rossouw et al. 2018). With these concepts in mind, many studies have committed to unravel the mechanisms that determine the dynamics of carbohydrate transport in plants (Münch 1927;Dewar 1993;Daudet et al. 2002;Thompson and Holbrook 2003a, b). Continuous monitoring of the carbohydrate flow within plants is impractical, making modelling a feasible approach to investigate the carbohydrate transport among various sources and sinks in plants ( De Schepper et al. 2013;Hall and Minchin 2013;Seleznyova and Hanan 2018).
The majority of carbohydrate allocation and transport models can be divided into two groups: (i) Common assimilate pool (CP) models, which assume that the ability of developing plant organs to attract carbohydrates is independent of the topological position of the organs (Brown et al. 2019). (ii) Mechanistic coupled phloem/xylem carbohydrate transport (CT) models, which assume that carbohydrate fluxes are proportional to the osmotically generated pressure differential at local scale, as induced by the activities of sinks and sources, e.g. carbohydrate unloading by different organs and rates of photosynthesis (Thompson and Holbrook 2003a;Allen et al. 2005;Seleznyova and Hanan 2018). The CP model has been widely used in various modelling platforms and is proven successful in simulating the biomass production and final yield of many annual crops (Heuvelink 1995;Brown et al. 2019). However, many studies involving tree or vine plant species show that organ growth on a particular branch is dependent on the carbohydrate status of that branch, in conjunction with the carbohydrate status of the whole plant (Piller et al. 1998;Pallas et al. 2010;Eltom et al. 2013). These results therefore suggest a semi-autonomy in terms of carbohydrate distribution, at branch scale (Pallas et al. 2016;Auzmendi and Hanan 2020).
To deal with the variability of shoot and fruit growth within the plant, transport models based on an electric circuit analogy have been developed for the peach tree (Allen et al. 2005) and the kiwifruit vine (Cieslak et al. 2011). At each developmental step of the system, the stationary phloem sap transport equations were explained using an electric circuit analogy with internodes represented as resistors and the carbohydrate sources and sinks associated with the lateral organs, represented by electromotive forces (Prusinkiewicz et al. 2007). However, the electric circuit analogy methods only consider the mass flux of carbohydrate rather than the overall flux of the phloem sap. As a result, the transport mechanism in these models cannot be interpreted as Münch flow, but rather a process similar to stationary diffusion (Seleznyova and Hanan 2018). To resolve this, Seleznyova and Hanan (2018) developed a coupled phloem/xylem transport model, which was the first transport model to provide continuous distribution of the system variables in a complex developing structure (Table 1). The model accounted for the non-linear dependence of phloem resistance and osmotic potential on the local carbohydrate concentration. The plant structure was modelled at a phytomer level with the Table 1. Summary of symbols used in the coupled phloem/xylem transport model.

Symbols
Definitions Unit c(x) Phloem carbon concentration at position x gC cm −3 Φ(x) Carbon potential defined as c 2 (x)/2 gC 2 cm −6 j Mass carbon transport gC h −1 R p A single (total) phloem resistance to phloem sap per unit length, estimated as a non-linear function of c(x) MPa cm −4 h R s Phloem resistance to carbohydrate mass flux per unit length of phloem vasculature gC h cm −7 r s Total internode resistance to the carbohydrate flux gC h cm −6 K max,i Maximum hydraulic conductivity of a stem segment gC s −1 m MPa −1 K i Actual hydraulic conductivity of a stem segment gC s −1 m MPa −1 ψ i Mean water potential of segment i MPa ∏ (c) Phloem osmotic potential MPa X p Internode 'xylem pull', a term representing the effect of xylem water potential difference on phloem transport gC 2 cm −6 K M,i Michaelis constant which defines the response of sink activity to Φ (x) and sink priority as well Unitless g (c) Source loading limiting function Unitless s i Structural carbon of organ i gC r Radius of a stem segment m l Length of a stem segment m ρ Density of a stem segment g m −3 M rsp,i Maintenance demand of organ i gC Comparison between the common assimilate pool and phloem carbohydrate transport models • 3 internodes represented by conduit elements and the lateral organs represented by sources and sinks. Transport equations were solved analytically for each internode before the solutions are adjusted and 'sewn' together using an iterative computational procedure, taking into account concentration-dependent sources and sinks. So far, the coupled phloem/xylem transport model has only been tested theoretically on an idealized sieve tube against analytical and numerical solutions obtained by Hall and Minchin (2013) and Thompson and Holbrook (2003a). To our knowledge, no studies have compared the results of phloem carbohydrate concentration (c(x), x represents the position in the transport pathway), carbohydrate allocation and organ dry matter (DM), as obtained by either the CP model or the CT model. This is partly because those two methods are fundamentally different and thus difficult to compare. The goals of the current study are to: (i) integrate the phloem/xylem transport model developed by Seleznyova and Hanan (2018) into the GrapevineXL model (Zhu et al. 2019); (ii) develop a statistical method for calibrating the CT model, thereby reducing the boundary for the adoption of CT models; (iii) calibrate the CT model against detailed experimental data and (iv) compare the performance of the CT model and the CP model under both homogeneous and heterogeneous canopy architectures by using the same carbon loading/unloading equations for both models (Fig. 1). A homogeneous canopy architecture is defined where all fruits were equally close to the carbohydrate sources or evenly distributed within the canopy, while a heterogeneous canopy architecture is defined as the distance between fruit and carbohydrate sources varies from one fruit to another fruit.

Model development
The models applied in the current study were based on GrapevineXL (Zhu et al. 2018(Zhu et al. , 2019, which is a GroIMP (Kniemeyer 2008)-based 3D functional-structural whole-plant grapevine (Vitis vinifera) model. GrapevineXL simulates the effects on post-véraison grape berry growth under a static canopy architecture, of variability in hourly environmental conditions (including radiation and soil water potential), plant water status (including the xylem and leaf water potentials) and plant carbon status. However, this model assumes uniformity in terms of the xylem water potential and phloem carbohydrate concentration c(x) (unit: gC cm −3 ) throughout the stems, without considering the hydraulic properties of the stem segments (internode, cordon and trunk) and the changes of c(x) along the carbohydrate transport pathway. To better understand the effects of canopy structure, and the distribution of c(x) along the stems on fruit growth and fruit variability, we made the following improvements: (i) clarifications of different carbohydrate sources and sinks based on the findings of Cieslak et al. (2011); (ii) adding the phloem/xylem transport model (Seleznyova and Hanan 2018); (iii) incorporation of the temperature response of all carbohydrate loading and unloading processes (Parent et al. 2010) and (iv) improving the canopy architecture representation. In this model, we use carbon atom as the unit for all carbohydrate calculations and mass balance checks, and the word carbon in following sections Figure 1. Illustration of the similarities and differences between the common assimilate pool model and a coupled phloem/xylem transport model. The same functions that define the carbon loading and unloading and their repose to c(x) for each organ are used for both models. For the common assimilate pool model, c(x) was calculated based on the assumption that carbon loading from leaves, stem and roots was equal to carbon unloading by stem, roots and berries at each hour. No distance and c(x) gradient along the transport pathway were considered. Here, the stem was just a simplified notation for all internodes (current season shoot), cordons (2-year-old shoot) and trunk (perennial woody part), although these objects were treated individually in the 3D model. For the coupled phloem/xylem transport model, c(x) was calculated at the bottom and top of each organ. Zhu et al. refers to the elemental carbon in the form of carbohydrate. The current model uses a static architecture as an example, but the methods are applicable to dynamic growth as shown by our preliminary tests.

Water transport.
Hydraulic properties of the stem segments were added based on functions described by Albasha et al. (2019). The maximum hydraulic conductivity of a single stem segment (K max,i , kg s −1 m MPa −1 ) was estimated based on the segment diameter (Tyree and Zimmermann 2013;Albasha et al. 2019). The actual stem hydraulic conductivity (K i ) varies with water potential because of the development of xylem embolisms under conditions of water deficit (Tyree and Zimmermann 2013).
where D i (m) is the average diameter of segment i, ψ i (MPa) is the mean of water potential of segment i. ψ * stem is the critical stem water potential that defines the steepness of the reduction in K i due to cavitation, and k 1 , k 2 and k 3 are dimensionless shape parameters.
The hydraulic conductance of the root and leaf was included in GrapevineXL and was estimated as a function of transpiration, circadian rhythms and [ABA] xyl . For a detailed description refer to Tardieu et al. (2015) and Zhu et al. (2018). The effect of xylem embolism on leaf hydraulic conductance was included as well (Zhu et al. 2018).

Carbon transport model.
The coupled phloem/xylem transport approach, developed in L-Studio (Seleznyova and Hanan 2018), was integrated into GrapevineXL to simulate carbon transport and within-vine competition (illustrated in Fig. 1 and Table 1). For solving the transport equations analytically, the coupled phloem/xylem transport approach introduced a new term called 'carbon potential Φ(x) (g 2 cm −6 )', which was defined as c 2 (x)/2. This new term is useful for expressing the mass carbon transport (j, g h −1 ) in a linear form (Seleznyova and Hanan 2018).
where R p (MPa cm −4 h) is the phloem resistance per unit length, which is estimated as a non-linear function of c(x), while R s (g h cm −7 ) is the phloem resistance to carbon mass flux per unit length of phloem vascular system. dψ(x)/dx is the xylem water potential gradient at position x. d ∏ (c)/dc (MPa g −1 cm 3 ) is the derivative of phloem osmotic potential to carbon concentration. R 0 (MPa cm −4 h) is the phloem resistance at zero phloem carbon concentration. R (0.0243 MPa cm 3 g −1 K −1 ) is the universal gas constant. T K is the absolute temperature. r 1 , r 2 , r 3 and r 4 are fitted parameters for calculating the increase of phloem resistance as a function of phloem carbon concentration. p 1 and p 2 are fitted parameters for calculating the increase of phloem osmotic potential as a function of phloem carbon concentration. Note that the current parameter for phloem osmotic potential and transport resistance were estimated based on the form of sucrose carbon, which is the most common carbohydrate form for transport in grapevine (Zhang et al. 2006). However, the mathematical and computational methods do not depend on the type of the main carbohydrate in the phloem and thus can be used for other carbohydrates as well, e.g. sorbitol. The carbon potential distribution in a segment with a length of l can be calculated using the following equations: where r s (g h cm −6 ) (equals to R s times l) is the total internode resistance to the carbon flux. X p is the internode 'xylem pull' representing the effect of the xylem water potential difference on phloem transport (subscripts b and t correspond to the bottom and the top of the internode).
Following the linearization of carbon potential distribution, the source/sink functions were linearized as well. Sink activities (carbon unloading from the phloem) were represented by a Michaelis-Menten function ( f (K M,i , Φ )) and then expressed in a linear form via the first two terms of the Taylor series.
where V max is the maximum carbon demand based on the organ's intrinsic properties and environmental factors (e.g. temperature). K M,i is the Michaelis constant, which is also useful to compare the priority in carbon unloading among sinks. The effects of K M on the rate of carbon unloading were shown in Supporting Information- Fig. S1. a and b are coefficients for the linearized form. Source activities, represented by carbon loading into the phloem, were estimated according to the potential available carbon for loading Comparison between the common assimilate pool and phloem carbohydrate transport models • 5 times a source loading limiting function ( g(c)), based on carbon potential in the vicinity.
where c 0 is the loading threshold (the inflection point) and g k is the response rate parameter. The effects of c 0 and g k on the rate of carbon unloading were shown in Supporting Information- Fig. S2. The value of c 0 was determined based on the average phloem sugar concentration for active loaders (0.21 %wt/wt in terms of sugar which is 0.098 %wt/wt in terms of carbon atom unit) as described by Jensen et al. (2013), and then converted into carbon potential (4e-3 g 2 cm −6 ). During simulation, at each time interval the model finds steadystate solutions for Φ (x) and j within the plant. GroIMP has builtin numerical methods for solving ordinary differential equations (Hemmerling et al. 2013), which in principle could be applied for solving mechanistic Münch-based equations for coupled phloem/ xylem transport in systems with dynamic structure; however, these methods were implemented mainly for diffusion-driven transport. Explaining Münch-based equations would add additional implementation complexity due to non-linearity of the problem. For instance, osmotic potential and transport resistance have non-linear relationships with phloem carbohydrate concentration. Hence, a special procedure was developed in GroIMP to force it to perform the same calculation sequences as per the L-Studio version (Seleznyova and Hanan 2018). We compared our results with the original model in the L-Studio under a single sieve tube with the same parameter and set-up [see Supporting Information- Fig. S3]. After verifying the results, we adapted the method on a real plant situation entailing numerous branching points.
The adapted procedure first linearized the sink/source activity of all organs (noted as branch flux), and then summarized carbon loading/unloading on a per-node basis starting from the end of each shoot and moving towards the point of attachment to the cordon/trunk at the base of the shoot. Subsequently, the same procedure was applied from the end of cordon/trunk to the root. It solved the carbon potential at the intersection between structural root (diameter > 2 mm) and fine root (diameter ≤ 2 mm) based on the summed fluxes. Afterwards, the carbon potential in each phytomer from root to shoot tip was updated in sequence, based on the carbon potential at the bottom of the phytomer and the total branch flux through the phytomer. Finally, a new carbon flux for each phytomer was calculated based on the local carbon potential and its own branch flux. This process was repeated until the summed error between the new carbon flux and carbon flux in the previous iteration for each phytomer was smaller than 1e −4 . Mass balance was checked at each time step to ensure the correct implementation of the method.

Clarifications of carbon sources and sinks.
The overall carbon demand was divided into the requirements for (i) structural growth, (ii) maintenance and (iii) the replenishment of reserves. The structural growth demand was further divided into demands for primary growth and secondary growth. Primary growth mainly refers to axial organ growth, while secondary growth mainly refers to the thickening of the organ, e.g. radial growth of the internodes, trunk and structural roots. Starch and total soluble sugars were combined to express the abundance of total NSC. The carbohydrate reserve synthesis and hydrolysis were calculated based on NSC. Since, in the current study, we use a static architecture of a post-véraison vine, only the primary growth of fine root and berry were considered; therefore, leaf and internode primary growth were not included in the simulation.
The root module as formerly used was divided into a fine root and structural root module because of their distinguishable functionalities (de Herralde et al. 2010). Fine roots are largely responsible for plant water and nutrient uptake and have a short lifetime. The fine root module has functions to calculate root primary growth, root turnover, carbon reserve abundance and hydrolysis, root length, soil to root surface resistance under different soil water potential, the xylem sap abscisic acid (ABA) concentration and root conductance. The structural root module represents the thicker roots, which are responsible for the anchoring of the plant in the soil and for carbon reserve storage, and have a longer lifetime. The structural root module has functions for calculating secondary growth (root thickening), root turnover, and carbon reserve abundance and hydrolysis.

Primary growth demand.
The primary growth of fine roots was simulated based on the method proposed by Cieslak et al. (2011), where its growth was relative to the structural biomass, and responded to atmospheric or soil temperature.
where s froot is the structural carbon of the fine root, k froot is the maximum relative growth rate of fine root, which is set to 1.25e-4 gC gC −1 h −1 based on the peak growth of kiwifruit (Actinidia deliciosa) roots during summer (Buwalda 1993). q fRoot g (0.2 gC gC −1 ) is the growth respiration coefficient. t froot (2e-5 gC gC −1 h −1 ) is rate of turnover of fine root (Buwalda 1993). The turnover rate for structural roots was assumed a quarter of that of fine roots (5e-6 gC gC −1 h −1 ) (Klein and Hoch 2015). f (T) is the temperature response of the primary growth. Parameters were determined based on the leaf expansion rate of kiwifruit vines under different temperatures measured in controlled growth chambers (Seleznyova and Greer 2001). The same temperature response was applied to root turnover, berry growth and secondary growth. A different temperature response was used for NSC synthesis and hydrolysis as these parameters were driven more by enzyme activities rather than cell division and elongation.
To simplify the calibration of carbon allocation, the primary growth of the berry was simulated using a logistic growth function (Eq. 16) instead of using a biophysical berry growth module such as originally implemented in GrapevineXL.
where the constant k berry defines the relative bunch growth rate and S berry,max is the potential bunch carbohydrate mass (in this model we simulate the bunch as whole instead of individual berries). s berry refers to the total bunch carbohydrate mass, that is, the mean berry carbohydrate mass multiplied with the number of berries.

Secondary growth demand.
The secondary growth of internodes, the trunk and structural roots was simulated using a constant relative growth rate multiplied with a temperature response (Cieslak et al. 2011). The rapid initial radius growth of the internode was excluded from the current model because we started the simulation at the post-véraison phase of grapevine development, after the conclusion of rapid vegetative growth. The carbon demand was then calculated by the radius (r, unit m), radius change (dr/dt, m h −1 ), length (l, m) and wood density (ρ, g m −3 ) (Cieslak et al. 2011).
where s i is the structural carbon of organ i. k sec is the long-term radial growth rate. k sec was first estimated based on increases in trunk circumference over a 14-year period, as determined for 800 field-grown Sauvignon blanc vines in New Zealand, and then optimized for the current experiment based on the change of structural root, trunk and shoot DM. The current model assumes the structural root as one long root, and its wood density as being the same as the density of the trunk. The length of the structural root was set as the length of the trunk multiplied by the ratio between the root structural biomass and the trunk structural biomass. As a result, the biomass ratio will stay relatively constant despite the effects of root turnover on the evolution of root structural biomass.

Maintenance demand.
Maintenance demand ( M rsp,i ) was modelled as the influx of carbon into a maintenance sink (Cieslak et al. 2011), with the responses to temperature and the ratio of NSC to structural carbon (SC) of the organ in question incorporated (Noguchi 2005).
where m i is the maintenance coefficient for organ i. m base is the minimum respiration percentage when the NSC to SC ratio is zero. Q 10 is defined as the increase in the respiration rate resulting from a temperature increase of 10 °C. The mean value of Q 10 (1.7) for leaves, bunches and stems, as measured by Poni et al. (2006) on 4-year-old potted Sangiovese grapevines, was used.

Reserve dynamics.
Carbon reserves as present in internodes, trunk, cordon and structural roots were modelled as active competing sinks driving starch synthesis. The rate of NSC synthesis depends on organ size and is limited by overloading. Remobilization of stored carbon is proportional to the amount of NSC in the organ and the rate of starch hydrolysis. The dynamics of carbon reserves are expressed by the following equation: where s res,i is the abundance of the reserve. k syn is the rate of NSC synthesis and k hyd rate of NSC hydrolysis. k syn and k hyd were constant for internodes, the trunk and structural roots, but different values were used for the leaves. For leaves, a fixed portion of carbon assimilates are stored as transitory NSC (12.5 %) and the rest are made available for loading into phloem during the day (Chew et al. 2014). A fixed portion (5 %) of the total NSC at each step is hydrolysed into sugars, available for loading during both the day and night, which eases the transition in leaf loading between the day and night. In the current simulation, we assume one structural carbon can hold one non-structural carbon based on the non-structural carbon concentration reviewed in Holzapfel et al. (2010).
is the ratio of enzyme activity at T k and at 293 K. f * (T k ) is the modified Eyring's equation proposed by Johnson et al. (1942). The numerator is the Eyring equation. The denominator (reversible denaturation of enzymes) is determined by enthalpy (∆HD) and entropy (∆SD) between the catalytically active and inactive states of the enzyme or enzymatic system. The ratio ∆H D /∆S D determines the temperature at which half of the enzymes are in the inactive state, and affects the temperature at which the rate begins to decline (Parent et al. 2010;Gauthier et al. 2020).

Canopy representation.
The canopy representation of GrapevineXL was improved in two major aspects: leaf angle and leaf shape.

Leaf angle.
Leaf angle includes leaf azimuth, that is, the leaf 's midrib angle in relation to the horizon and roll angles around the midrib. The leaf angle changes in reaction to local radiation conditions and its distribution would vary greatly under different training systems. To account for leaf angle, we adopted the concept of turtle optimization method, developed by G. Louarn (INRAE, Lusignan, France, pers. comm.). The leaf angle optimization procedure started from the highest leaf rank on each shoot progressing downwards for each leaf.
Comparison between the common assimilate pool and phloem carbohydrate transport models • 7 During the optimization, each leaf was replaced by a semi-hemisphere (named as 'turtle') consisting of 72 radiation sensors (Supporting Information-Video S1). The radiation sensor does not interfere with radiation transport, but only measures the surface irradiance (Kniemeyer 2008). The semi-hemisphere was positioned according to the global axis and put on a horizontal plane. The centre of the semihemisphere was positioned at the leaf insertion points on the shoot, while the radius of the semi-hemisphere was set to the length of the petiole. Afterwards, the radiation model was evoked and the virtual radiation sensors recorded the radiation intensity from each direction. The corresponding leaf was then rotated according to the direction of the radiation sensor that exhibits the maximum irradiance among the different sensors. This approach accounts for leaf position and orientation within complex grapevine canopies.
2.1.4.2 Leaf shape. The leaf surface was changed from 2D surface to 3D surface based on the Virtual Riesling model (Schmidt and Kahlen 2018;Schmidt et al. 2019). Furthermore, as the grapevine leaves are relatively large (~100 cm 2 ) and there is considerable variability in radiation intensity within one leaf, the 3D leaf surface was further divided into 16 triangles, termed leaf-facets. Photosynthesis and transpiration were calculated in each small leaf-facet based on the local radiation intensity and then summarized to estimate the whole of leaf photosynthesis and transpiration rates.

Carbon allocation data set
The carbon allocation between the different grapevine components was calibrated based on a detailed potted vine study involving three levels of leaf availability, during the berry ripening period. The leaf availability treatments corresponded with 100 leaves retained per vine, 25 leaves retained per vine and no leaves (i.e. all leaves removed at the start of the experiment) (Rossouw et al. 2017a). The study was started 9 days after the beginning of berry ripening, i.e. 9 days after the onset of véraison (i.e. berry softening) ( Fig. 2; see Supporting Information- Fig. S4). Forty own-rooted cv. Shiraz (clone EVOVS12) grapevines, grown in 30-L pots, containing commercial potting mix, were used in the study. The grapevines were enclosed in an outdoor bird-proof cage in the hot climate Riverina grape-growing region in New South Wales, Australia and were well-watered throughout the growing season. The 3-year-old grapevines were spur-pruned to five two-bud spurs in the winter and distributed in four rows with 10 vines each. Three vines from each treatment were destructively harvested every 9-10 days after the initiation of the treatments. For each grapevine, the total leaf area, fresh and DM of the whole-root system, leaf blades, trunk, stems, leaves and all berries were determined. The total NSC and nitrogen concentrations were determined for the roots, leaves and berries. For NSC analysis, soluble sugars were first extracted from a 20 mg subsample of each tissue using 3 × 1 mL × 10 min washes of 80 % aqueous ethanol. The first two volumes were at 80 °C and the third at room temperature ). After centrifuging between each wash, the three aliquots were combined, diluted to 10 mL and the concentration of sucrose, d-fructose and d-glucose determined with commercial enzyme assays (Megazyme International, Bray, Ireland). For starch analysis, the remaining wood sample was resuspended in 200 µL dimethylsulfoxide and heated at 98 °C for 10 min. The remainder of the analysis was then performed using commercial enzymes and glucose assay kits (Megazyme International). Briefly, 300 µL thermostable α-amylase in MOPS buffer was added, mixed and incubated for 15 min in a 98 °C water bath. After cooling, 400 µL amyloglucosidase in sodium acetate buffer was added and incubated at 50 °C for 60 min. The samples were mixed at 20-min intervals, and then centrifuged at 10 000 rpm for 2 min. Supernatant from root samples was diluted 1:11, and leaf samples 1:6 in Ultra-pure water. Glucose concentration of the diluted samples was then determined colorimetrically and the amount of starch in the original 20 mg sample calculated.
The mass balance at the start of the experiment was investigated in the no leaves treatment by comparing the initial abundance of NSC in all organs and DM increases (mainly in fruit) after 38 days. Unfortunately there was a considerable gap between those two values that can be explained by the possible contribution of DM from amino acids or protein etc. The gap could arise from the loss of fine and structural root components during the destructive harvesting process, underestimation of the NSC concentrations in the fine and structural root by the sampling and measurement methods, and variations in organ DM from vine to vine. Such a discrepancy between DM increases and carbohydrate analyses has been documented in the past. Jordan et al. (2000) found the discrepancy between total soluble fraction measured by chemical analysis and by refractometer could not be accounted for from typical soluble acid concentrations reported in the literature, or by other identified sources of error. Further investigations are required to understand the underestimation of carbohydrate concentration and content by the sampling and analysing methods. As the carbon allocation method requires mass balance as a prerequisite for usage of the data set, the fine root and structural root DM and NSC concentrations were systematically increased by 30 % for the whole data set and for all treatment. The NSC concentrations in the trunk and shoot stem were not measured during the trial, and were set to 50 % of the corresponding values in structural roots, as based on the measurements in a similar pot trial (Rossouw et al. 2017b). This roughly made the initial NSC reserves in all organs equal to the final biomass at the end of experiment in the no leaves treatment. Although the above-mentioned factors may affect the relative values of the parameters we obtained, the comparison of the CP and CT models shall not be affected since the same set of parameters are used across the two models.

Simulation set-up and model initialization
The GrapevineXL model was revised to have a similar canopy architecture compared to the vines in the leaf availability trial ( Fig. 2; see Supporting Information- Fig. S3). The virtual vine had five short spurs attached to the trunk. Each spur bore two shoots and each shoot bore one bunch. The weight of each bunch was calculated as the observed mean bunch weight per vine at the initiation of the leaf availability treatments, divided by 10. The leaf size profile of the Cabernet Sauvignon vines, which was parameterized in GrapevineXL, was scaled to produce a similar leaf area as measured in the leaf availability trial, described above (Zhu et al. 2018). Potential berry growth rates, initial DM and NSC concentrations in each organ, leaf area and the dynamics of the leaf nitrogen content were determined based on the results of the leaf availability trial, described above (Rossouw et al. 2017a). The rate of photosynthesis was determined by a CO 2 response curve measured on a similar potted vine experiment, based on the same variety and environmental conditions. Air temperature, radiation, relative humidity and wind speed were obtained from the Wagga Wagga Amo station (35°09′29″S, 147°27′27″E), which is about 14 km from the experimental site. Soil water potential was set to −0.02 MPa to ensure no water stress was present. Leaf angles were optimized at the start of the each simulation.
Parameters related to carbohydrate unloading, e.g. V max and Michaelis constant, were first taken from the GrapevineXL (Zhu et al. 2019) and L-Kiwi (Cieslak et al. 2011) models. The Michaelis constants were further scaled to match the conversion of carbon concentration to carbon potential in the equation, and then explored by trial and error to ensure that the simulated trends were in agreement with experimental data. Parameters related to phloem osmotic potential and resistance were derived from Seleznyova and Hanan (2018).
Simulations were made with one plant but cloned nine times to remove border effects, through the GridClonerNode function in GroIMP. The nine plants were configured in three rows with three plants in each row. Row and plant distance were both set to 1 m. Radiation absorption by each leaf was calculated through a GPU-based raytracing method provided by the GroIMP platform (Henke and Buck-Sorlin 2018). For simplicity, only diffused radiation was used to represent the light environment. Diffused radiation was estimated using an array of 72 directional surface light sources positioned regularly in a hemisphere of six circles with 12 light sources each (Chelle and Andrieu 1998;Xu et al. 2011;Zhu et al. 2015). The intensity of the hourly total radiation was interpreted from the recorded daily total radiation assuming a diurnal sine shape.

Sensitivity analysis
Eight key parameters that determine carbon allocation between the different grapevine organs were selected for global sensitivity analysis and optimization (Fig. 3). These parameters are the rate of NSC hydrolysis ( k hyd ) and synthesis (k syn ), the secondary growth rate (ksec), the Michaelis constant for carbohydrate unloading towards NSC synthesis (KM,res), secondary growth (KM,sec), berry growth (K M,berry ), fine root growth (KM,froot) and the response rate for the source loading function (gk). Given the high number of parameter-based combinations and the computational time needed for one simulation (10 min × 3 leaf treatments: 100 leaves retained per vine, 25 leaves retained per vine and no leaves), Figure 2. Illustration of the model configuration for different leaf number treatments corresponding to the pot experiment by Rossouw et al. (2017a) [see Supporting Information- Fig. S3], and the effects of leaf angle optimization on leaf position (panels A-C). An animation of how leaf angle optimization was done from top rank to low rank was shown in Supporting Information-Video S1. The potted vine was spur-pruned to five two-bud spurs in the winter and trained to have 10 shoots per vine. During simulation, the observed mean crop load at the start of the leaf treatment was equally divided into 10 bunches and distributed to each shoot on the third internode.
Comparison between the common assimilate pool and phloem carbohydrate transport models • 9 a fractional factorial design was used to reduce the number of parameter combinations (Lecarpentier et al. 2019). Three levels were considered for each parameter. A fractional factorial design of resolution V (3 5 = 243 simulations) was chosen to ensure estimation with no confusion of main effects and pairwise interaction of input factors. Third-order and higherorder interactions were excluded. The fractional factorial design was generated via the planor package in R software (Monod et al. 2012).
The generated parameter set was loaded into the GrapevineXL model and the model simulated for all combinations of the parameter set and leaf treatments (243 × 3 simulations). Each simulation took about half an hour on a Plant & Food Research computer cluster node with 32 cores. The simulation results were then analysed to determine the sensitivity of the selected model outputs to each parameter, through the R multisensi package (Bidot et al. 2018). Principal components analysis with two axes was used for dimensional reduction. Main sensitivity index (MSI) and interaction sensitivity index (ISI) were included as outputs for each parameter.
where a and b represent different parameters, SS a is the sum of squares associated with the main effect of a. SS T is the total sum of squares of the output. SS a,b is the interaction sum of square for parameter a and b.

Parameter optimization
The simulation results (243 * 3 simulations) were further used for developing Gaussian process emulators, that is, a statistical model linking the parameter values and the model output as fitted using the R package GPfit (MacDonald et al. 2015). An emulator was developed for each combination of leaf treatments, selected model outputs and days after leaf treatments when the measurements were conducted. One observed value requires one emulator. The emulators were used to optimize the parameter values in accordance to the observed data set. The optimization was undertaken separately on all combinations in terms of the leaf availability treatments, through the R DEoptim package (Ardia et al. 2011). To balance the contribution of the different treatments to the parameter value, mean values across different combinations of the treatments were used. The criteria for optimization were set as minimizing the sum of squares between observed values and predicted values.
((X sim,j,i − X obs,j,i )/X obs,j,i ) 2 (26) Figure 3. Diagram for the global sensitivity analysis and parameter optimization. Firstly, 243 parameter sets were generated using the R planar package with three levels for each parameter and resolution V (3 5 = 243). All parameter sets were entered into the GrapevineXL model and the model was run for all combinations of parameter sets and leaf treatments (243 × 3 simulations).
The simulation results were then analysed to determine the sensitivity of the selected model outputs to each parameter through the R multisensi package. The simulation results were further used to develop statistical emulators using the R package GPfit. An emulator was set up for each combination of leaf availability treatments, selected model outputs and days after leaf treatment initiation, i.e. when the measurements were done. One emulator corresponds to one observed value. The emulators were used to optimize the parameter values given the observed data set. The optimization was done through the R DEoptim package.
where i is the sample number or time point, n is the total number of measurements, j is the variable number, m is the total number of variables applicable, X sim,j,i is the simulated value for variable j and sample number i and X obs,j,i is the observed value for variable j and sample number i. Bunch DM, trunk DM and total root NSC (fine root NSC + structural root NSC) were used for optimization. Due to the large variability in total root DM, it was excluded for optimization. Sensitivity analysis and parameter optimization were performed based on the CT model; the protocols are illustrated in Fig. 3.

Scenario simulations
In contrast to the homogeneous canopy architecture as applied for the modelling as relevant to the leaf availability trial, a heterogeneous canopy architecture was used for further understanding the effects of the proximity of bunches to the carbohydrate sources, in addition to the crop load, on berry growth. Therefore, a single cane-pruned system was used, but with only one shoot with 10 leaves at the end of the cane. The crop load was varied by adding completely defoliated shoots to the cane, with each shoot bearing two bunches. In total, seven different crop load scenarios were simulated, ranging from 4 to 16 bunches per vine (noted as treatment 1-7). Initial plant conditions and weather conditions were the same as in the leaf availability trial, except for the architecture of the shoot. The bunch weight and c(x) at each position were used as the outputs of the simulation.

Statistical index
Goodness-of-fit between observed and simulated values was assessed by root mean square error (RMSE): where i is the sample number, n the total number of measurements, X sim,i is the ith simulated value and X obs,i is the ith observed value. The units of RMSE are equal to those of the observed data. Standard deviation (σ, Eq. 25) and coefficient of variation (CV, Eq. 26) were used for checking the variations in the simulated bunch DM under heterogeneous canopy scenarios.
where µ is the treatment mean of the simulated values.

Sensitivity analysis
Sensitivity analysis allowed for the variance of the DM of bunches, trunks, structural roots and structural root NSC to be separated into main effects, second-order interactions and residuals, which represent the higher-order interactions and uncounted variances ( Fig. 4; see Supporting Information- Fig. S5). Total main sensitivity indices for the eight parameters were above 0.95 for each output under all leaf availability treatments. Interaction sensitivity indices were relatively small compared with MSI (<5 %) and residuals were all smaller than 0.011. The sensitivity of the model output to certain parameters changed under different leaf number. For the 100 leaves per vine treatment, where the carbohydrate supply was relatively abundant, the bunch DM was most sensitive to the secondary growth rate (ksec, MSI 0.44), followed by K M,berry (0.25) and K M,sec (0.16). When no leaves were present, and the carbohydrate supply was therefore severely restricted, bunch DM was most sensitive to the rate of NSC hydrolysis (k hyd , 0.41), followed by K M,res (0.18), K M,berry (0.15), k sec (0.10) and K M,sec (0.09). For the 25 leaves per vine treatment, the contributions of different parameters to the total variance were, however, more equally spread.
Structural root NSC was most sensitive to k hyd , and its sensitivity increased steadily with a reduction in the leaf number per vine, ranging from 0.38 to 0.81 from the 100 leaves per vine to no leaves treatment. In contrast, the sensitivity of k syn decreased from 0.14 to 0.02 with a reduction in leaf number per vine. K M,res played an important role in regulating the structural root NSC with an average MSI of 0.15 across the three treatments. Structural root NSC was very sensitive to k sec (0.29) under the 100 leaves per vine treatment; however, it became insensitive under both the 25 leaves per vine treatment and when no leaves were present. A similar pattern was found for structural root DM. However, a contrasting pattern was found for the trunk DM [see Supporting Information- Fig. S5], which was most sensitive to K M,res (mean MSI 0.39), followed by k hyd (0.21) and k syn (0.18).
In addition, the model output was sensitive to the construction of canopy architecture. Leaf angle optimization increased the whole-plant radiation interception by 14 % and carbon assimilation by 6 % across the 38-day simulation period [see Supporting Information- Fig. S6]. Dividing a whole leaf into 16 small facets decreased the overall carbon assimilation by 1 % as compared with the undivided leaf [see Supporting Information- Fig. S7]. The remaining results were obtained on simulations with leaf angle optimization and leaf-facet method.

Model calibration
The eight parameters in the global sensitivity analysis were optimized based on the workflow described in Fig. 3. After the optimization, g k, K M,froot and k syn were fixed to their optimized value, since these parameters are relative insensitive and k syn is correlated with k hyd . The remaining five parameters were used to generate a new 243 parameter sets which was subjected to a second round of optimization. The distribution of the model outputs of the new 243 × 3 simulations covered most of the DM and NSC variations as observed in the different leaf availability treatments [see Supporting Information- Fig. S8].
Large variations in the optimized parameter values were found for the different treatments and additionally for combinations of the treatments [see Supporting Information- Table S1]. The CV for the different parameters ranged from 0.26 to 1.37. The mean optimized values across the different treatment combinations were entered in the GPfit emulator, which successfully captured the variations in bunch DM (RMSE = 37.1, R 2 = 0.94), total root NSC (RMSE = 6.3, R 2 = 0.8) and trunk DM (RMSE = 6.8, R 2 = 0.52; see Supporting Information- Fig. S9).

Model verification and carbon allocation
The mean optimized parameters were entered in GrapevineXL to verify the results. The model was able to capture the dynamics of berry DM, total root NSC and trunk DM for the 25 leaves per vine and no leaves per vine treatments (Fig. 5). However, it underestimated the total root NSC and overestimated the trunk DM for the 100 leaves per vine treatment. The model reproduced the dynamics of total root DM for the 25 leaves per vine treatment. It underestimated the total root DM under 100 leaves per vine, probably because of a low initial DM input and big variations in root DM between the first and second sampling date (Fig. 5).
Berries represented the largest carbohydrate sink among the different grapevine organs. In fact, approximately 80 % of the available carbon was allocated to the berries at the start of the experiment, irrespective of treatment (Fig. 6). However, the exact rate of carbon unloading into the berries ranged from 2500 mg C day −1 (sum of hourly carbon flux) for the no leaves treatment to 7500 mg C day −1 for the 100 leaves per vine treatment [see Supporting Information- Fig. S10]. The percentage of carbon unloaded by the berries decreased over time for the 100 leaves per vine treatment. This is largely attributed to the decrease in carbon demand from the berries over time (Fig.  5), while the carbon demand from the stem, structural roots and fine roots stayed relatively stable (Fig. 6). The relative decrease in terms of the percentage of carbon unloaded by berries was much slower for the 25 leaves per vine and no leaves per vine treatments, compared with the 100 leaves per vine treatment (>70 % for the 25 leaves and no leaves per vine treatments compared with 25 % for the 100 leaves per vine treatment, at the end of simulation; Fig. 6).
Leaves contributed to more than 90 % of the total carbohydrate source for the 100 leaves per vine treatment, and for more than 50 % for the 25 leaves per vine treatment (Fig. 6). The percentage of carbon loading by structural roots increased from ~4 % for the 100 leaves per vine treatment, to ~25 % and ~57 % for the 25 leaves per vine and no leaves per vine treatments, respectively. No differences were found in the percentages of carbon unloaded or loaded by the different grapevine organs for both the CP and CT model, irrespective of the intact number of leaves per grapevine.

The dynamics of phloem carbohydrate concentration under homogeneous architecture
As expected, diurnal fluctuations were found for c(x) with maximum values found around the middle of the day, whereas minimal values were found around sunrise [see Supporting Information- Fig.  S11]. The diurnal pattern was strongly influenced by the air temperature [see Supporting Information- Fig. S11] and by the rate of leaf photosynthesis and carbon loading. c(x) decreased rapidly shortly after the implementation of the treatments (from the start of the simulation), particularly for the 25 leaves per vine and no leaves per vine treatments, coinciding with a reduction in air temperature (Fig. 7). Overall c(x) increased over time for the 100 leaves per vine treatment and decreased slowly for the no leaves per vine treatment (Fig. 7). It stayed relatively stable for the 25 leaves per vine treatment, regardless of the fluctuations caused by changes in air temperature.
Only a small gradient was found in c(x) from the apical parts of the plant towards the fine root for the 100 leaves per vine treatment, and no gradient was found for the 25 leaves per vine treatment (Fig. 7). The gradient of c(x) was, however, reversed for the no leaves per vine treatment, with positions close to the fine roots and trunk exhibiting the largest c(x) (Fig. 7). The c(x) values obtained by the common assimilate pool model were generally found to be within or slightly below the gradient for all three of the leaf availability treatments.

The effects of heterogeneous architecture on berry growth and c(x)
Large variations in bunch DM and c(x) at different positions of the vine were illustrated by the CT model under heterogeneous canopy conditions (Figs 8 and 9). The CV of bunch DM increased from 0.01 under treatment 1 (4 bunches per vine) to 0.17 under treatment 7 (16 bunches per vine). The DM of bunches on the leaf-bearing shoot (B1 and B2) decreased slowly with an increase in crop load (from 59 to 45 g), whereas the DM of other bunches decreased more sharply to a mean DM of 31 g under treatment 7 (Fig. 8). Interestingly under treatment 7, bunch 15 had the largest DM among all bunches except bunches 1 and 2, probably because it was the bunch closest to the trunk and structural root. The uniform bunch DM obtained by the CP model was statistically equal to the mean bunch DM value obtained by the CT model for each treatment (Fig. 8).
The mean c(x) at the middle of Day 3 decreased with an increase in the number of bunches per vine, while the variation of c(x) increased under the same scenario (Fig. 9, panel A). The daily maximum c(x) difference among all positions increased rapidly at first and then stabilized in conjunction with an increase in the number of bunches per vine (Fig. 9). These values ranged from 0.025 to 0.057 g cm −3 , which was close to the daily mean c(x) value found for the 25 leaves per vine treatment (Fig. 7).

Common assimilate pool versus coupled phloem/xylem transport
We compared the performance of CP and CT models using c(x) and organ DM as quantitative measures under different canopy architectures. Under a homogeneous canopy architecture, where the predominant carbohydrate sinks during the berry ripening period, that is, the different bunches, are all located the same distance from the carbohydrate sources. For this homogeneous canopy-based scenario, the gradient of c(x) along the carbon transport pathway was relatively small ( Fig. 7; see Supporting Information- Fig. S11). This result is in agreement with the findings of Heuvelink (1995), where a common assimilate pool was found to be suitable in terms of simulating tomato truss weight. The total amount of carbohydrate in the plant was found to be the only factor determining final tomato truss weight in their study. Similarly De  showed that the rate of carbohydrate loading per day, as estimated by a leaf level-based photosynthesis model, was close to the rate of carbohydrate loading as shown by a mechanistic flow and storage model in tomato plants. The Comparison between the common assimilate pool and phloem carbohydrate transport models • 13 flow and storage model related variations in stem diameter, measured at three different parts of the plant, with phloem carbohydrate loading and carbohydrate concentration dynamics in tomato. Under a heterogeneous canopy architecture, however, the gradient of c(x) and variation of bunch DM increased greatly with increasing crop load. This result demonstrates the benefits of using a CT model, instead of a CP model, for capturing the variations in c(x) under conditions of carbon source limitation, and particularly when variability exists in regards to the distance between carbon sources and sinks. In the study of Heuvelink (1995), as they pointed out, the plants were probably sink-limited. Under conditions where the sink size and/or strength is limited, the growth rate of all organs would be expected to be close to their potential rates, irrespective of the effects of c(x) on growth. Pallas et al. (2010) conducted six experiments on two grapevine varieties to quantify the effect of variations in carbohydrate supply and topological distances between sources and sinks on organogenesis, morphogenesis and biomass growth. They found that the CP model was inadequate for describing grapevine development, and that dividing the whole plant as a sum of independent axes could be a possible way of simulating biomass partitioning (Pallas et al. 2010;Auzmendi and Hanan 2020). Future studies could apply the model used in the current study to the experimental conditions as outlined in Pallas et al. (2010), in order to quantify the extent of c(x) variation in the system, and to help define the best approach for simulating the phenotypic plasticity. Figure 6. Simulated daily mean percentage of carbon loading into phloem by leaf, fine root, stem and structural root under different leaf treatments (panels A, C and E) and percentage of carbon unloading by berry, fine root, stem and structural root (panels B, D and F). Stem was a simplified notation for all internodes (current season shoot), cordons (2-year-old shoot) and trunk (perennial woody part). Results were obtained from the coupled phloem/xylem transport (CT), and no significant difference was found on the percentage of carbon loading and unloading by the grouped organ types between the common assimilate pool (CP) model and CT model in this simulation. Figure 8. The effects of heterogeneous architecture and crop load on bunch weight. The setups of the simulation for treatments 1 and 7 were illustrated in panel A, and the resulted bunch dry matter for each treatment was shown in panel B. In total, seven different crop loads (T1-T7) were simulated with two bunch increments for each treatment, ranging from 4 to 16 bunches per vine. Bunches were counted from left (on the leafy shoot) towards the trunk. Simulations were conducted under both the common assimilate pool (CP) and the coupled phloem/xylem transport (CT) models. Coefficient of variation was only calculated for the CT model by dividing the standard deviation by the treatment mean. The presented bunch weight was the value at the end of simulation (38 days after leaf treatment).

Carbohydrate allocation interacts with phloem carbohydrate concentration and distances between sources and sinks
The hierarchy in carbohydrate allocation among various sinks has been observed in several studies (Wardlaw 1990;Pallas et al. 2010). We adopted the Michaelis-Menten method as previously used in L-kiwi (Cieslak et al. 2011) to capture the hierarchical behaviour. The model illustrated that the percentage of carbohydrates allocated to the berries increased under source-limited conditions, that is, for the 25 leaves per vine and no leaves treatments (Fig. 6). This result is attributed to the fact that the rate of carbohydrate unloading interacted with c(x), and a low value of K M for berry (Table 2) suggests that the rate of reduction in carbohydrate unloading was slower with decreased c(x) for berry compared with other organs. In our model, the following hierarchy was configured: maintenance respiration > berry growth > fine root growth > long-term secondary growth = reserve synthesis, which is consistent with the findings of Pallas et al. (2010) and Rossouw et al. (2017b). Primary growth, e.g. leaf and internode elongation, was not quantified in the current study; however, primary growth is expected to exhibit a higher priority than berry growth (Cieslak et al. 2011).
The global sensitivity analysis showed that the DM of different organs was more sensitive to the sink activities, i.e. the parameters controlling the rate of carbohydrate unloading, for the 100 leaves per vine treatment, where the assumption is that no source limitations existed ( Fig. 4; see Supporting Information- Fig. S5). While under sourcelimited conditions, that is, the 25 or no leaves per vine treatments, the organ DM was more sensitive to source activities, particularly NSC hydrolysis. The comparison of the sensitivity index of trunk DM and structural root DM, which were both subjected to the same NSC hydrolysis and synthesis, and model parameters, suggests that the sensitivity is dependent on initial conditions. This suggests that the model could be used to identify the most limiting factors for carbohydrate accumulation and distribution, and may therefore help improve productivity in an authentic plant system.
The effects of the distance between the sources and the sinks on carbohydrate unloading were clearly demonstrated by the scenario simulation (Figs 8 and 9). Interestingly, the effect of distance proved much more pronounced under a higher crop load, particularly where variability in the bunch DM occurred (Fig. 8). This is partly explained by the increase in the gradient of c(x) along the pathway under a high crop load (Fig. 9), and partly because bunches 1 and 2 were directly on the leaf-bearing shoot and all other bunches were on separate branches (Fig. 8). An extra simulation was conducted to confirm the second hypothesis (i.e. the relative location of bunches 1 and 2) by moving bunches 1 and 2 to a separated shoot, close to the leaf-bearing shoot instead of directly on the leaf-bearing shoot. This reduced the advantage of bunches 1 and 2 in terms of carbohydrate unloading and final DM, and their mean weight was consequently reduced from 45 to 38 g. These results indicated the interaction between phloem carbohydrate concentration and fruit distribution in affecting the carbohydrate allocation among individual organs. Such an interaction cannot be readily captured by the common assimilate pool model even with the effect of source-sink distances on carbon partitioning (Pallas et al. 2016). The model developed by Hall and Minchin (2013) can provide analytical solutions for steady-state coupled phloem/xylem flow. However, their model needed to be developed for a particular architecture and cannot be applied to dynamic architecture. In contrast, the current method was readily applicable to the dynamics of growing perennial fruit crops. Figure 9. The effects of heterogeneous architecture and crop load on the mean and variation of c(x) within the pathway during the middle of Day 3, after the start of the simulation (A), and the daily maximum c(x) difference across the 38 simulation days (B). In total, seven different crop loads were simulated with two bunch increments for each treatment, ranging from 4 to 16 bunches per vine. The variation of c(x) within the pathway was assessed by exporting all c(x) on the main pathway, that is, c(x) associated with internodes, cordons, trunks, and structural and fine roots. The maximum hourly c(x) difference within the vine was calculated for each day and its daily variation was plotted in panel B. Blue dots in panel A represent the c(x) values as simulated by the CP model, and the red points were the mean of the CT model in both panels. The black line in the box plot represents the 50th percentile. The bottom and top of the box represents the 25th and 75th percentiles, respectively. The line above the box represents the values within 1.5 times interquartile (75th percentile to 25th percentile) above the 75th percentile. The line below the box represents the values within 1.5 times interquartile below the 25th percentile. Black points were the points not within the above defined ranges.

Future perspectives
Temperature showed a strong positive effect on c(x), especially for the no leaves treatment [see Supporting Information- Fig. S11]. This was largely related to the temperature response of enzyme activity, e.g. the enzymes involved in NSC hydrolysis and synthesis, which increases with temperature (Parent et al. 2010). However, Field et al. (2009) andField et al. (2020) showed that low soil temperature could stimulate starch accumulation in grapevine roots, which was demonstrated to be related Parameters were estimated in four complementary methods: (i) directly estimated from experimental data described above (experiment); (ii) directly taken from literature; (iii) taken from literature first but then adapted for grapevine based on the trends published in literature or in our data collection (exploration); (iv) taken from literature first but then calibrated for our data through numerical optimization (calibration).
Comparison between the common assimilate pool and phloem carbohydrate transport models • 17 to the dormancy response and root-synthesized cytokinins. Further work is required to quantify the full effects of air and soil temperature on carbohydrate dynamics, especially around the dormancy period. Water stress has pronounced effects on the connectivity between the xylem and phloem. Water stress changes both sink/source activities and thus can either increase or decrease the c(x). For example, water stress reduces the demand for primary growth, including leaf development, and also reduces photosynthesis. Furthermore, water stress increases the xylem water potential and reduces the water flow from the xylem to the phloem, thus increasing the density of the phloem sap and increasing c(x) (Savage et al. 2016). The current model can capture the effects of water stress on xylem water potential and photosynthesis (Zhu et al. 2018(Zhu et al. , 2019. However, the assumption of zero water flux resistance between xylem and phloem and water equilibrium between xylem and phloem may not hold under water stress conditions. Furthermore, the current carbon transport model focuses more on the c(x) gradient, which is the driver of carbon flux, rather than the absolute value of c(x).
The current study did not have detailed photosynthesis measurements to quantify the rate of carbon assimilation at the whole-plant level. This added a source of error in estimating the amount of carbohydrate supply under the 100 and 25 leaves per vine treatments. We did compare the simulated rate of transpiration per hour with the loss of pot weight to ensure our simulations were in good range. The loss of pot weight was measured in 2 h after irrigation around the middle of the day, approximately 106 ± 19 g h −1 for no leaves per vine, 223 ± 16 g h −1 for 25 leaves per vine and 461 ± 6 g h −1 in a sunny day. More dedicated studies where we can monitor the rate of carbohydrate assimilation, the types of the carbohydrates and their corresponding concentrations in the phloem vessel adjacent to different organs would help us better understand the carbohydrate metabolism in plants and their responses to environmental conditions. Incorporation of leaf angle optimization increased the total carbon assimilation by 6 % under our outdoor pot environment, while adopting the leaf-facet method decreased the carbon assimilation by 1 % compared to whole-leaf method [see Supporting Information-Figs S6 and S7]. The reduction caused by the leaf-facet method was expected as the rate of increase in photosynthesis (μmol m −2 s −1 ) decreases with increasing light intensity (Zhu et al. 2018). When dividing a single leaf into small facets, some facets capture additional light to the detriment of other facets. However, the net positive gain of the illumined facets will not outweigh the net loss from non-saturated facets. Further work could evaluate the effects of the leaf-facet method and leaf angle optimization under more heterogeneous light environments with sun-flecks and partial leaf shading and more realistic canopy architecture and training systems, e.g. vertical shoot positioning and Geneva double curtain training systems.

CON CLUSIONS
The current study compared the common assimilate pool model with the coupled phloem/xylem transport model within a whole-plant functional-structural grapevine model. The results show that under a homogeneous canopy architecture, where the fruit was evenly distributed, the performances of the two models were very similar. However, under a heterogeneous canopy architecture, where the distance to the main carbohydrate sources varies for each fruit, noticeable differences between the two models are found. The coupled phloem/xylem transport model offers greater potential in terms of understanding the variation in c(x) and fruit DM, which is a goal of many horticultural practices in vineyard and orchard. Furthermore, our whole-plant model showed that the most limiting factor for fruit growth is dependent on the plant source/sink status, and the model can, therefore, help with disentangling the effects of different processes on fruit growth, and offer practical suggestions for vineyard and orchard management.

SUPPORTIN G INFOR M ATION
The following additional information is available in the online version of this article- Table S1. The optimized parameter values for carbon allocation based on different combinations of leaf treatments. Video S1. Leaf angle optimization animation. Figure S1. The effect of K M values on the rate of sink unloading corresponding to carbon potential and carbon concentration. Figure S2. The effect of inflection point of the source loading function c 0 and scaling parameter g k on rate of source loading. Figure S3. Comparison of the rates of carbon flux and carbon potential outputted by the carbon transport procedure implemented in GroIMP and the original model in L-Studio under a single 30-cm sieve tube. The sieve tube was divided into 300 nodes with 1 cm for each internode. Figure S4. Illustration of the potted vines that used in the experiment of Rossouw et al. (2017a). Figure S5. Global sensitivity analysis of structural root dry matter and try dry matter to eight parameters under 100 leaves per vine, 25 leaves per vine and no leaves treatment. Figure S6. The effects of leaf angle optimization method on the total light absorption and carbon assimilation over the 38 simulated days. Figure S7. The effects of leaf-facet method on the total light absorption and carbon assimilation over the 38 simulated days. Figure S8. Distribution of the simulated bunch dry matter (DM) and total root non-structural carbon, total root DM, trunk DM in the 243 simulations for parameter optimization in the second round. Figure S9. Verification of the simulated bunch dry matter, total root non-structural carbon, trunk dry matter by the GPfit emulator. Figure S10. Simulated daily rates of carbon loading into phloem by leaf, fine root, stem and structural root under different leaf treatments and carbon unloading by berry, fine root, stem and structural root. Figure S11. The diurnal dynamics of c(x) under different leaf treatments.

SO URCE S OF FUNDIN G
This work was part of the Plant and Food Research Wine Research programme, funded by the Ministry of Business, Innovation and Employment (MBIE) Strategic Science Investment Fund. The work of M.H. was supported by the European Regional Development Fund-Project 'SINGING PLANT' (no. CZ.02.1.01/0.0/0.0/16_026/0008446).