Coupling kinetic models and advection-diffusion equations. 1. Framework development and application to sucrose translocation and metabolism in sugarcane

The sugarcane stalk, besides being the main structural component of the plant, is also the major storage organ for carbohydrates. Previous studies have modelled the sucrose accumulation pathway in the internodal storage parenchyma of sugarcane using kinetic models cast as systems of ordinary differential equations. To address the shortcomings of these models, which did not include subcellular compartmentation or spatial information, the present study extends the original models within an advection-diffusion-reaction framework, requiring the use of partial differential equations to model sucrose metabolism coupled to phloem translocation. We propose a kinetic model of a coupled reaction network where species can be involved in chemical reactions and/or be transported over long distances in a fluid medium by advection or diffusion. Darcy’s law is used to model fluid flow and allows a simplified, phenomenological approach to be applied to translocation in the phloem. Similarly, generic reversible Hill equations are used to model biochemical reaction rates. Numerical solutions to this formulation are demonstrated with time-course analysis of a simplified model of sucrose accumulation. The model shows sucrose accumulation in the vacuoles of stalk parenchyma cells, and is moreover able to demonstrate the up-regulation of photosynthesis in response to a change in sink demand. The model presented is able to capture the spatio-temporal evolution of the system from a set of initial conditions by combining phloem flow, diffusion, transport of metabolites between compartments and biochemical enzyme-catalysed reactions in a rigorous, quantitative framework that can form the basis for future modelling and experimental design.


Introduction
Sugarcane is part of the family of grasses (Moore and Maretzki, 1996). Unlike some of the other grass crops, like maize, sorghum and barley, it is the stalk (also called a culm) which acts as the primary sink for carbohydrates produced during photosynthesis (Komor, 2000b). The stalk is divided into alternating nodes and internodes, with a single leaf attached to a node. Sucrose is primarily synthesised in the leaves and loaded into the phloem, where it is transported up or down the stalk and unloaded in the storage parenchyma (Komor, 2000a;van Bel, 1993;Walsh et al., 2005). Sucrose is the most abundant soluble carbohydrate found in sugarcane and is actively accumulated to levels much higher than in most other plants (Bull and Glasziou, 1963; Moore and Maretzki, 1996).
Kinetic modelling is a useful tool for analysing the control and regulation of metabolic pathways (Rohwer, 2012), and the application of systems biological approaches in the context of sugarcane physiology have been reviewed in Rohwer and Uys (2014). To understand better the factors affecting and controlling sucrose accumulation in the sugarcane culm storage parenchyma, Rohwer and Botha (2001) have proposed a kinetic model of the most important biochemical reactions involved. A key result from this model is the important role played by the futile cycling of sucrose (Whittaker and Botha, 1997), i.e. the process by which sucrose is broken down and resynthesised, which is known to expend metabolic energy in the hexokinase and fructokinase reactions. The model showed that factors leading to a decrease in futile cycling cause a concomitant increase in sucrose accumulation. A variation of the model was subsequently used by Schafer et al. (2004) to study the kinetics various isoforms of the enzyme sucrose synthase. Bosch (2005) added the branch towards trehalose synthesis to the original model and concluded that partitioning to trehalose does not significantly impact sucrose accumulation. Uys et al. (2007) extended the original model of Rohwer and Botha by explicitly modelling the sucrose synthase and fructokinase isoforms, as well as adding the phosphofructokinase, pyrophosphate-dependent phosphofructokinase and UDP-glucose dehydrogenase steps. In an attempt to model physiological changes during culm maturation, a steady state was calculated for eight different sets of enzyme maximal activity data, representing internodes 3-10. Importantly, all the above models only focused on the reactions around sucrose and only in the cytoplasmic space (symplast) to study the effect of sucrose breakdown and synthesis; they did not explicitly include other subcellular compartments.
Overall, the model in Uys et al. (2007) showed higher sucrose concentrations than that in Rohwer and Botha (2001). However, these levels still did not come close to measured concentrations, most probably because the vacuole, apoplast and phloem were not explicitly included. To address these shortcomings and develop a more realistic model of sucrose accumulation, we present in this paper a framework for combining phloem flow and reaction networks into a single model. Explicitly modelling the phloem compartment requires an entirely new approach, because now geometry and sap flow become important and the reactions can no longer be assumed to occur in a "well-stirred" homogeneous space. The processes involved in sucrose accumulation can be collectively referred to as advection-diffusion-reaction (ADR) behaviour, and as shown below, this can be modelled as a set of partial differential equations (PDEs).
Since the ADR framework uses PDEs, classical methods for determining and quantifying the effect of a parameter change on system behaviour (such as metabolic control analysis; Kacser and Burns, 1973;Heinrich and Rapoport, 1974) can no longer be applied. To overcome this, and to provide two solutions for performing a parameter-sensitivity analysis on such systems, the Morris sampling method and the Fourier Amplitude Sensitivity Test (FAST) are introduced in the accompanying paper (Uys et al., 2021).

Modelling advection-diffusion-reaction systems
We first review briefly the existing approaches for modelling phloem flow, and then integrate them with biochemical reactions to formulate the ADR framework.

Models of phloem flow
Sap flow in plants involves both the xylem and the phloem and has been studied extensively for many decades; a number of models have been formulated and a recent comprehensive review (Jensen et al., 2016) focuses on the biophysical aspects. The phloem transports carbohydrates (mainly sucrose) from the sites of photosynthesis to the sink tissues in the plant, and comprises a whole complex of cells. The sieve elements are arrayed end-to-end, separated from each other by a sieve plate, forming the sieve tube. The sieve tube in turn is surrounded by companion cells, one or more to each sieve element (van Bel, 2003;van Bel et al., 2002). Thompson (2006) described phloem with three analogies.

"A sieve tube is like dialysis tubing"
There is constant movement of solutes and water between the phloem and the surrounding tissue. Turgor pressure is regulated through changes in solute concentrations. If phloem loading and unloading were only to occur at the extremities of a phloem tube, pressure gradients would be difficult to control as local changes in pressure gradients become difficult to induce.

"A sieve tube is like a capillary"
Phloem flow is dominated by shear forces rather than inertial forces; in other words, viscosity trumps momentum. As a consequence, pressure gradients develop as a result of energy dissipation. If viscosity increased (i.e., viscous forces started to dominate), the flow rate of solute would decrease. In this context, Thompson and Holbrook (2004) studied "information transmission" through the phloem. Under certain conditions it is possible for a local change in solute concentration or pressure to propagate through the phloem faster than the actual velocity of the phloem sap (Thompson, 2005(Thompson, , 2006. The implication is that message passing along the phloem is possible in a reasonable amount of time.

"A sieve tube is like a full bathtub"
A change in the solute concentration at a local point along the phloem translates into a flux of solute, towards or away, from the surrounding sap. This ripple effect propagates along the phloem tube and becomes more diffuse as it moves further away from the source of the perturbation. The absolute magnitude of the flux is proportional to the solute concentration and also the gradient of the solute concentration. Bieleski (2000) points out four features of sugarcane that separate it from other plants. 1. The plant can be seen as a row of sinks all connected by the phloem.
2. Sinks with widely differing sugar content all "feed" off the same phloem strands.
3. Phloem unloading occurs along the entire length of the phloem tube. This is in contrast to other plants where the phloem has a well defined "terminus" in the roots. 4. Any point along the phloem can act as a loading or unloading site.
The MÃ¼nch hypothesis (M• unch, 1926, 1927, 1930) (Knoblauch and Peters, 2017, provide an historical overview) is the premise for most modelling of fluid flow in xylem and phloem (Daudet et al., 2002;Lampinen and Noponen, 2003;Holtta et al., 2006). Fluid flow is driven by a pressure gradient generated by an osmotic potential; the effect has been confirmed by Gould et al. (2005). Thompson and Holbrook (2003a) have developed a mass-balanced non-steady-state model of phloem flow and in Appendix A of their paper also provide a concise review of phloem modelling based on Münch's hypothesis. De Schepper et al. (2013) have reviewed mechanisms and controls of phloem transport, and more recently, Minchin and Lacointe (2017) have formulated a model to analyse the effects of phloem unloading and reloading on flows between source and sink.
In a series of follow-up articles, Holbrook (2003b, 2004) used dimensional analysis to study their approach to phloem modelling. The assumption of equilibrium between phloem sap and apoplast water potential was shown to be plausible (Thompson and Holbrook, 2003b). This is important for the model formulation presented later, since passive water flux between the phloem and apoplast need not be accounted for (Thompson, 2006). Furthermore, a relatively small pressure drop over a small distance would allow more accurate control of solute loading.

Transport phenomena
ADR models are examples of a class of problems called transport phenomena (Beek et al., 1999). The mathematics describing the movement of a measurable quantity turns out to be similar for a large range of phenomena. Central to any transport phenomenon is the principle of conservation. For example, the diffusion of heat or the diffusion of a solute through a liquid would both depend on the divergence of the gradient together with some, possibly constant, coefficient. Mathematically the expressions are identical even though two different physical properties are being modelled. The differences are the actual values and units of the variables and parameters, together with the boundary and initial conditions imposed. All of these can be measured experimentally.

The ADR framework
For a summary of the notation used, refer to Table 1 S is involved in a set of enzyme catalysed reactions,  , found in  , with each reaction having signed stoichiometric coefficients, i c , and rate, v .

2.
S crosses O because of (a) convection, e.g. advection by a velocity vector field, (b) diffusion due to a concentration gradient, with some coefficient, D , or (c) facilitated or active transport due to a set of transporters,  , with each transport step having signed stoichiometric coefficients, the following set of PDEs desribes a system in which a species, ij s , may be advected, diffused, transformed into another type of species or be moved into a neighbouring compartment, where u is the phase average velocity and e  is a unit vector.
Define the del operator ) , , ( where the general rate vector, To solve Equation (3) for s the following need to be specified: the equations that govern reaction rates and phase average velocities, diffusion coefficients, the geometry and subsequent domain on which the variables are defined, the initial conditions and possible boundary conditions that may exist. Given sufficient initial conditions the time-dependent behaviour of the system can be modelled.
As it stands, Equation (3) is a very broad framework in which a model needs to fit. Not all compartments will allow advection to occur, some species may diffuse quickly through a compartment to form a homogeneous concentration field, other compartments will allow some combination of advection, diffusion and reaction to occur. The next section deals specifically with the simplifications that can be made to Equation (3) to model translocation of photoassimilate in the phloem.

"Solute driven advection"
Long-distance photoassimilate transport in plants occurs via the phloem , Bel(2003) Error! Reference source not found.. All further modelling in this work will only consider a single solute (i.e. sucrose) in the phloem, therefore we will only consider Equation (3) for one variable. The right hand side of Equation (3) can also be reduced to phloem loading and unloading steps, which yields ).
The phase average velocity, u  , can be obtained from Darcy's law , where  is the hydraulic conductivity (a measurable quantity assumed to be constant,]ThoHol03b,  is the viscosity, P is pressure and P  is the pressure gradient. By substituting Equation (6) into the second term of Equation (5) the following is obtained, This formulation casts advection in terms of a pressure gradient, which can in turn be solved by using the van't Hoff equation for dilute solutions and de~Paula(2002) Error! Reference source not found., where  is the osmotic pressure and 0 P is the surrounding pressure from the apoplast.
Substitution of Equation (8) into Equation (7) gives the following relation, The conservation equation is therefore of the following form, This allows phloem translocation to be modelled with the following approximation: ). ( called this "solute driven advection".
The next section describes the equations governing catalytic conversions and crossmembrane transport, in other words, finding the rates to populate v in terms of metabolite concentrations and parameters.

Reactions and cross-membrane transport
where s , p and m are respectively substrate, product and modifier concentrations scaled by their specific half saturation constants. f V is the limiting forward reaction rate,  is the mass action ratio and eq K is the equilibrium constant for the reaction (the disequilibrium ratio eq /K  determines the direction of the reaction); h is the Hill coefficient and Cornish-Bowden(1997) Error! Reference source not found., ;  and  are the scaling coefficients, either activating or inhibiting the enzyme through, respectively, a K-effect or V-effect . The dot product, r m , is used for convenience.
The equation for uni-uni reactions is still simpler: Downloaded from https://academic.oup.com/insilicoplants/advance-article/doi/10.1093/insilicoplants/diab013/6206762 by guest on 02 April 2021 For the dynamic system as defined in Equation (3), enzymes are also subject to movement by advection-diffusion. This means that enzyme concentrations (call them e ) will change with position and therefore their limiting rates (maximal activities) will also change. Nothing prohibits treating the enzyme concentrations as state variables as well with appropriate equations for enzyme synthesis and degradation in v . However, in this work it is assumed that the enzyme concentration gradient remains constant with time, i.e.

Compartmental modelling
When solutes are transported from one compartment to another the difference in compartment volume needs to be accounted for. A possible way to do this is to have volumes appear explicitly in the relevant equations, with solutes measured in amounts. Terms in the differential equations would then have units of amounts per unit time. Therefore all s in Equations (12)- (14) are actually defined as (15) and appear as such in the model. Solute flow and reactions are defined on certain compartments and sometimes also span a number of compartments.

Justification for phenomenological equations
Modelling sometimes requires a choice to be made between using equations of state that are mechanistically complete or just behaviourally complete, the latter being a necessary condition for the former. In general, the definitive equation for modelling fluid flow would be the Navier-Stokes equations et~al.(1999)Beek, Muttzall, and van Heuven Error! Reference source not found.. There are many difficulties in solving these equations and what is gained from the solution is not necessarily of any practical use. For example, the Navier-Stokes equations can account for turbulent flow, but one would not expect to see turbulence in plant phloem. In other words, Reynolds numbers would be very low and flow would be laminar. Similarly, phloem sap is not really compressible. This would justify the use of simplified forms of Navier-Stokes as in et~al.(2002)Henton, Greaves, Piller, and Minchin Error! Reference source not found.. Similarly, Equation (11) is used as a simplified version of the conservation equation for modelling phloem .
A similar reasoning applies to enzyme kinetic rate equations. Mechanistic rate equations often require a greater number of parameters than phenomenological ones; moreover, these parameters have not necessarily been measured. The generalised reversible Hill equation introduced above explicitly incorporates a thermodynamic term and thus implicitly satisfies the Haldane relationship, has all kinetic terms separated from thermodynamic terms, and possible allosteric terms are separated from mass action terms. Finally, all kinetic parameters are phenomenological and thus have a direct operational definition in terms of their experimental interpretation For these reasons we have taken a purely phenomenological approach in the formulation presented above, i.e. the Darcy equation for modelling flow in porous media to model phloem flow, and the generalised reversible Hill equation to calculate reaction rates.

Computational methods
D is a generic coefficient of diffusion (for example heat or solute diffusion) and  S is a source term (i.e. a source of  ), then there is a natural way to program the equation in Python . It is clear that Equation (3) is amenable to this formalism.
Simulation data were visualised with the matplotlib plotting library . The aim of this paper is to combine phloem modelling including loading and unloading along the entire length of the sieve tube with compartmentation and detailed reaction kinetic modelling.

Pathway
The model is constructed to emulate a single sugarcane stalk. A segment of the whole pathway, a single node/internode pair, is shown in Figure 2. The primary storage organ is the stalk (in particular the vacuoles of the parenchymal cells) and sources of assimilate are distributed along the stalk. There are five nodes and internodes, which represent immature to mature tissue. The primary storage organelle is the vacuole, with possible assimilate remobilisation. A futile cycle known to exist in sugarcane is modelled explicitly; sucrose can be broken down into fructose and glucose. The hexoses can then enter the hexose phosphate pool and sucrose can be synthesised again. However, the cycle can either take place in the symplast alone (i.e. the cytosolic space inside the plasma membrane formed by linking of the cytoplasm of successive cells through plasmodesmata), or across the symplast and vacuole. The apoplast (i.e. the extracellular space outside the plasma membrane formed by cell walls of adjacent cells) was not included as a separate compartment in this first version of the model. The model also includes two moiety-conserved cycles ( is not meant to refer to protons in a mechanistic way, but is included to generically couple the uptake of S into the vacuole to a secondary concentration gradient. In our opinion this level of abstraction is justified in the core model presented here, the purpose of which is primarily to demonstrate the feasibility of the modelling approach. Metabolic activity can regenerate the moieties (e.g. reaction 11 in Figure 2 found., and Botha(1997) Error! Reference source not found., where accumulation is due to energy driven anti-porter proton gradients and Komor(1985) Error! Reference source not found., et~al.(2006)Reinders, Sivitz, Hsi, Grof, Perroux, and Ward Error! Reference source not found.. Stored assimilate can also be remobilised from the vacuole . Energy is captured from metabolite breakdown and consumed by synthesis and transport steps, reactions in the respiration branch are subject to negative allosteric feedback (i.e. synonymous with PFK and PK), and these reactions are positively cooperative.
The inclusion of moiety-conserved cycles such as ATP/ADP in a model can sometimes lead to the introduction of constraints on the fluxes because the moiety must be balanced at steady state. Such constraints may be unintentional or artefactual, and one way of counteracting a possible sideeffect is to introduce an additional ATP consuming or producing reaction. This was the reason for including reaction 12. Table 2 summarises the processes that were accounted for in each compartment of the combined model. Advection-diffusion is assumed only to occur in the phloem and is dominated by the advection component, as derived in Equation (11) above. Species in the symplast can diffuse, be transported or react; source and vacuole by contrast are assumed to be "well-stirred". Thus, depending on which compartment is being modelled, the appropriate terms from Equation (3) for the various processes are included or ignored.

Partial differential equations
Keeping the above in mind, the equations governing the system are divided into four groups. The first group describes those species that are only involved in reactions, in other words species found in the leaf and vacuole: the PDEs governing species behaviour in the symplast are: Note that the stoichiometry for reaction 10 is The stoichiometric model in Figure 2 has two conservation relations describing the total amounts of The i v appearing on the right hand side of the PDEs are themselves functions and are determined by the rate equations (Equation (12)).

Geometry and domain
The first step in building the model was to describe the geometry, which needed to be divided into finite volumes or be "meshed". Compartments were then defined on these elements. A compartment can span a number of elements and compartments of the same type can be separated by a number of elements. A pathway or reaction network was defined across the compartments. Species concentrations and reaction rates are variables that in turn were defined on the network. Note that the same chemical species found in two different compartments was modelled as two separate variables. Therefore each variable was only defined for its specific compartment and could not assume a value anywhere outside of it (the domain of the variable). In the model, reaction rates were calculated from species concentrations and it follows that they cannot have values at positions where a species variable does not exist. The behaviour of the variables is governed by a set of PDEs. Once the model formulation had proceeded this far, a time course simulation could be performed given sufficient parameter values, initial conditions and boundary conditions. The model behaviour was then analysed with appropriate software and programs.
Given the above hierarchy of modelling steps, we started with a straight line segment to approximate the domain by a 1d mesh with arbitrary length, , and divided into 50 finite volumes. The compartments, i.e. source (leaf), phloem, symplast and vacuole, were then defined on these finite volumes (Figure 3). Each finite volume element has at least 3 compartments, i.e. the phloem, symplast and vacuole. Some finite volume elements have an additional fourth compartment, i.e. the source. The phloem and symplast span all 50 elements, with leaves acting as source tissue feeding assimilate into the phloem at the nodes, and storage parenchyma in the internodes (symplast) acting as sink tissues extracting assimilate from the phloem (see Figure 3). There are five sequential node and internode pairs, numbered from 1 to 5, where internode 1 is considered immature tissue and internode 5 is mature tissue. There are ten elements to a node/internode pair; the first element of each group of ten is used for the node, the remaining nine for the internode. The only path out of a leaf element is into the phloem. There are four possible paths for a molecule in the phloem, up, down, into the leaf or into the symplast. Molecules in the symplast can move up, down, into the phloem, or into a vacuole. A vacuole corresponds to one finite volume; consequently there are 50 vacuolar compartments in the symplast. Molecules in a vacuole cannot move to another vacuole, they can only move to the symplast.
The choice of 50 = N was determined by the minimum number of finite volumes required to have five node and internode pairs such that the length of an internode was 9 times longer than a node, as approximation of real internode lengths Dillewijn(1952) Error! Reference source not found. .
Because this was a core model, any units assigned to parameters and variables would be arbitrary to some extent. However, it should be emphasised that variables were defined in amounts, so concentrations-and by implication also half-saturation constants-were in amount per volume.
Rates take their units from the maximal activities and these were in amount per unit time, not concentration per unit time. Each of the finite volumes was assumed to have identical size across the different compartments, obviating the need for correction to kinetics for movement between volumes of different size. For the compartments allowing advection and/or diffusion (phloem and symplast), the simulated concentration of a species assumed the same value throughout a particular finite volume; this is the consequence of the discretisation of the continuous PDEs.
The "ground state" for the model had all parameters set to a value of 1-so-called "vanilla" kinetics. Any changes to the ground state were only introduced to account for a certain behaviour. Furthermore, parameters were chosen to exaggerate the desired behaviour, as discussed below.

Thermodynamic parameters
Equilibrium constants were chosen with the following points in mind: 1. the preferred direction in which a reaction was to occur, 2. adherence to the Haldane relationship, 3. adherence to microscopic reversibility and detailed balance, and 4. dependencies between equilibrium constants for coupled reactions. Point 1 was addressed by not assuming any preference for substrates or products, except in the case where "energy equivalents" were produced or consumed. The use of generic reversible Hill equations implicitly satisfied Point 2.
The concepts of detailed balance and microscopic reversibility require that the product of the equilibrium constants around a closed cycle be equal to 1. Points 3 and 4 are closely related. Consider The easiest way to satisfy all these criteria was to specify the equilibrium constants for the half reactions as follows:

Kinetic parameters
All substrate and product half saturation constants were set to 1. For reactions 10 and 11, 0.8 = =   and the modifier half saturation constant was 5 . The Hill coefficient in both reactions was set to 4 , because both PFK and PK are tetramers. In reality the Hill coefficients are frequently lower, as the number of subunits represents an upper bound on the numerical value the Hill coefficient can take. However, as mentioned, the idea behind this model was to provide a simplified, if slightly exaggerated, representation with which the behaviour of the system could be simulated in a reasonable amount of CPU time.
Amongst the maximal activities in the model that remained constant with increasing z are those of assimilate synthesis in the leaf and of the phloem loading step. They were chosen to be the highest activities, because there are only 5 points along the stalk that produce assimilate through phloem loading, but 50 that consume it through phloem unloading. The chosen value was 10 = 1,2 f V .
The second group of activities chosen to remain constant were those of the two reactions involved in the futile cycle, the active proton transport across the tonoplast, and the remobilisation step from the vacuole-reactions 4, 5, 7 and 9. These were kept in the vanilla state, i.e. 1 = .
The maximal activities that were modelled to change with increasing z are shown in The rough guiding principle in selecting parameters was to observe accumulation of S along the stalk and with time. The ground state parameter set (all values 1 = ) unsurprisingly did not show this, and we tried to induce this behaviour with as few parameter changes as possible.

Flow coefficients
Diffusion coefficients were set at the same value of 5 10  for all species in the symplast. The flow coefficient was assigned a value of , 10 (20) such that the diffusion coefficient in the symplast was at least one order of magnitude smaller than the flow coefficient in the phloem. Since diffusion and translocation occur in opposite directions with reference to a concentration gradient, the convention is to choose negative flow parameters, so that fluxes are positive. In other words, if the concentration gradient has negative slope, then the direction of flow is positive. The initial values of each species variable, given in amounts, are summarised in Table 3. The amounts of the moiety-conserved species add up to a constant sum as follows:

Boundary and initial conditions
10,

Simulation results
This section summarises the development of concentration and reaction rate profiles during the time simulation. The results are presented in Figures 5 (concentration profiles) and 6 (reaction rates).

Sawtooth profiles
The distributions of P along the stalk all exhibited "sawtooth" profiles to some extent. These peaks correspond to the node positions; this is where all the assimilate is being loaded and phloem unloading rates are the highest. In the model, species are advected (in the phloem) or diffuse (in the symplast) away from the nodes. If the rate of translocation, advection or diffusion is slower than the reaction rates or rate of facilitated transport, then species form a heap round the nodes.

"Sucrose"
The highest concentration throughout the simulation was that of the "sucrose" analogue S , as can be seen on the left hand side of Figure 5.
S concentrations were initially evenly distributed along the stalk, and end the simulation at a shallow gradient monotonically decreasing with z . Similarly S in general decreased down the stalk, except for the local peaks around the nodes. The S concentration profile developed from an initial linearly increasing distribution to one that became constant with z -once gain with the exception of the node peaks. Towards the end of the simulation, the concentrations of S decreased from leaf, to phloem, to symplast and then increased again in the vacuole. In other words, assimilate flow was down a concentration gradient.
S increased along the stalk, with concentrations at immature end climbing faster than those at the mature end.

Moieties
The simulation profiles of moiety pairs,  v had a very pronounced negative parabolic profile that gradually disappeared within the first 150 time-points. This is discussed in detail in the Discussion section and Figure 8.

Futile cycling
The reactions 4 v and 5 v form a futile cycle in the symplast. The calculated rate of 5 v was higher than 4 v , therefore more S was being broken down than regenerated. Moreover, S is also broken down in the vacuole by 6 v , adding to the overall degree of futile cycling as illustrated in Figure 6. Note that in spite of this futile cycling ratio of 1 > it was still possible for S to increase over time, because S is additionally also taken up into the symplast through 3 v . We investigated whether this phenomenon could be captured with our simplified core model. To model the effect of shading, the following event was triggered at 102 = t during a timedependent simulation of the reference model.
In other words, all but the second leaf had the maximal velocity for reaction 1 set to zero (Figure 7). The following effects can be noted after the "shading" event, compared to the reference model: v (sucrose synthesis in the leaf) and 2 v (phloem loading) in the second (i.e. nonshaded) leaf increased above levels of the reference model. They remained higher for the rest of the simulation.
• P accumulated more slowly than in the reference model.
• Metabolite concentration profiles decreased more rapidly with z compared to the reference model.
• Assimilate flowed up and down the stalk away from the only remaining source.

Discussion
In this paper we have presented a general advection-diffusion-reaction framework for modelling phloem flow in plants coupled to biochemical reactions and transport across compartment boundaries. The framework was applied to a simplified model of sucrose accumulation in sugarcane. While the modelling of both phloem and xylem flows has a long history DisplayText cannot span more than one line!, and recently has been extended to a whole-plant framework ?, termed CPlantBox,]ZhoSchVan20 for modelling water and carbohydrate flows, to our knowledge this is the first report combining explicit simulation of enzyme-catalysed reactions within the ADR framework.
In the following we first discuss the mathematical framework, followed by a critical evaluation of the model description and further discussion of some simulation results.

Framework
Equation (3) is a general framework that allows the modelling of metabolites that are involved in chemical reactions, advected or diffused. It is a broad framework, which requires that an expression for fluid velocity be provided, compartmentation be defined, rate equations and stoichiometry be specified, and a suitable geometry with boundary conditions be provided. This was successfully applied to a pathway of biochemical reactions. By allowing such a pathway to vary spatially, more complicated behaviour can be modelled compared to a well-mixed case. Furthermore, the model provides a good first approximation of the modelling of sugarcane geometry.
There are, of course, a number of limitations in this approach. First, the use of Darcy's Law to model fluid flow could be interpreted as an oversimplification. Darcy's Law is a phenomenological equation that describes fluid flow in porous media. It requires flow parameters to be measured. The use of Darcy's Law is justified where macroscopic flow is modelled, i.e. it is entirely inappropriate for modelling fluid flow in a single sieve element, but adequate for modelling flow through many sieve elements strung together. A statistical averaging effect is at work over sufficiently large scales, and phloem flow and translocation over a whole plant qualifies as such. Despite its limitations, the use of Darcy's Law has advantages: for example, it is possible to simplify Equation (3) by taking into account the mechanism of phloem translocation . This simplification is easier to solve, because, first, the convective term is removed, and second, no pressure gradient has to be calculated. The simplification also describes the macroscopic (or phenomenological) aspects of phloem translocation, i.e. that flow is due to a pressure gradient induced by a solute.
The use of the van't Hoff equation is a further simplifying assumption. This equation is usually a first approximation for the calculation of osmotic pressure, since it only applies to dilute solutions. However, the only viable alternative is to use an interpolative function for a specific solute using known experimental measurements. The usual way of increasing the accuracy of the van't Hoff equation is by using a virial-like expansion for non-ideal cases, but the virial coefficients need to be measured by experiment and de~Paula (2002)  . and Minchin(2008) Error! Reference source not found. calculated that the error between a first and higher order approach is in the order of 8%. Such an error is insignificant for the simplified model with arbitrary parameter values presented here, and we therefore did not adopt the more complicated approaches.

Model
The model formulation in Figure 3 is of course a simplification of what happens in a sugarcane plant, and the formulation could be criticised on this ground. In what follows, we provide a justification for our simplifying choices. To model growth a method such as "adaptive meshing" would be required to dynamically adapt the domain to the change in geometry et~al.(2003)Plewa, Linde, and Weirs Error! Reference source not found.. PDEs that account for volume changes and partitioning of carbohydrates to structural components would be needed as well; however, this was beyond the scope of out initial "proof-of-concept" of the ADR framework.
Third, there is an achievable method of modelling growth with the current formulation, namely by altering the Neumann boundary conditions. Suppose . The modelling approach thus far has been to set the gradient of all state variables to zero at the boundaries, ensuring zero flux across the boundaries. Although changing the fixed gradient to a non-zero value would not explicitly model growth, the supply or demand for reactants could then take place across the boundaries. Furthermore, this exchange can change with growth. This is left for future work. Fourth, the model discretisation in terms of a 1d -mesh is of course a simplification. For the sugarcane stalk, which is unbranched, this is an appropriate approximation; the same argument could be made for constructing similar models for other grasses. The approximation is, however, insufficient for other plants that show branches in their stems. In this case, a 2d -or 3d -mesh would have to be defined that appropriately approximates the physical structure to be modelled. While this is tractable in principle, and the FiPy software has built in support for various 2d -and 3d -meshes as well as custom meshes generated by Gmsh and Remacle(2009) Error! Reference source not found., a custom 3d finite element mesh generator, this comes of course at a vastly increased computational cost. To a lesser extent the same applies to M and N . It is not as big an issue because both these species exist in the same compartment. Diffusive effects even out the concentrations to a uniform distribution if a local rate were to temporarily violate the moiety sum. Sawtooth profiles are a result of phloem loading being faster than translocation away from the source. Since assimilate flow is down a concentration gradient, this means that flow can be up or down the stalk, which is particularly evident for the sock experiment (see below).

Sock experiment
The model response of The model thus shows how "photosynthetic" rates may respond in the absence of genetic regulation. Further model adaptations could look at the combined effect of mass-action responses and genetic upregulation of an enzyme. For example, the accumulation of P is slower in the attenuated model compared to the reference model. An increase in the maximal velocity of reaction 1 (genetic upregulation) would allow the attenuated model to compete more closely with the reference model.
Finally, the metabolite and rate gradients as seen in Figure 7 may not reflect the situation in the plant completely realistically. The local accumulation of a metabolite reflects points to a diffusion or advection limitation; a higher diffusion coefficient would allow faster movement away from the source.

Conclusions
We have presented a general advection-diffusion-reaction framework for the integrated spatio-temporal simulation of plant physiological processes, which couples phloem flow to diffusion of species, biochemical reactions, and transport across compartment boundaries. The framework was applied to a simplified model of sucrose accumulation in sugarcane, which is able to capture the spatio-temporal evolution of the system from a set of initial conditions. Despite the simplicity of the model, a number of experimentally observed features of plant metabolism could be reproduced. The quantitative framework presented can form a rigorous basis for future modelling and experimental design, and also allows for the identification of key controlling steps through various forms of sensitivity analysis, as shown in the accompanying paper et~al.(2021)Uys, Hofmeyr, and Rohwer Error! Reference source not found..     The sugarcane stalk was approximated in 1d and compartments defined as indicated. Finite volumes of identical size are defined on each compartment; for the phloem and symplast these are shown by dotted lines because solutes can be transported along these continuous compartments by advection and/or diffusion as specified in Table 2.  S is only defined at a node and therefore appears as "bars". Those concentrations with a "saw-tooth" profile each have a small peak corresponding to a node position. Note that all the S species showed high concentrations. Note the parabolic profile of reaction 6, this is discussed in Figure 8 and the text.   A possible explanation for the prevalence of higher rates in the centre of the stalk (e.g. 6 v in the reference model, see Figure 6) is provided by increasing maximal activities with L upstream of the considered rate, and decreasing maximal activities with L downstream form it. See text.

Data
• adr_normal_data.npy Data for baseline model evaluation.

COMPETING INTERESTS
There is NO Competing Interest.