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

Abstract

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 upregulation 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.

1. 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 synthesized in the leaves and loaded into the phloem, where it is transported up or down the stalk and unloaded in the storage parenchyma (van Bel 1993; Komor 2000a; 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).

A number of factors influence the amount and rate of sucrose accumulation. Environmental conditions such as temperature, sunlight (or rather incident radiation), rainfall and soil types play a role (van Dillewijn 1952; Alexander 1973). Agricultural practices such as fertilization and chemical ripening can speed up accumulation (van Dillewijn 1952; Alexander 1973). Biochemical factors such as genetics (Ming et al. 2001), futile cycling of sucrose (Wendler et al. 1990; Komor 1994; Whittaker and Botha 1997; Zhu et al. 1997; Rohwer and Botha 2001) and membrane transporters will limit the maximum accumulation possible (Lemoine 2000; Rae et al. 2005b; Reinders et al. 2006).

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 resynthesized, 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 Schäfer 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 uridine diphosphate (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).

2. METHODS

2.1 Modelling ADR systems

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

2.1.1 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 et al. 2002; van Bel 2003). 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, 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ünch 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; Hölttä 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, Thompson and 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.

2.1.2 Transport phenomena.

Advection–diffusion–reaction 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.

2.1.3 The ADR framework.

For a summary of the notation used, refer to Table 1. Suppose a molecule, $S$⁠, with concentration, $s$⁠, is found in a small control volume, $Θ$⁠, bounded by a surface, $O$ (Fig. 1). The number of $S$ in $Θ$ can change if the following occurs:

1. $S$ is involved in a set of enzyme-catalysed reactions, $ρ$⁠, found in $Θ$⁠, with each reaction having signed stoichiometric coefficients, $ci$⁠, 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, $ci$⁠, and a rate, $v$⁠, as in the case of a reaction.

Let $i∈σ$ index a set of species, $j∈χ$ index a set of compartments, $k∈ρ∪τ$ index a set of reactions and transporters, and $l∈{1,2,3}$ index a set of Cartesian coordinate axes. Using the following shorthand,

the following set of PDEs desribes a system in which a species, $sij$⁠, 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 $∇≡e→1∂x+e→2∂y+e→3∂z=(∂∂x,∂∂y,∂∂z)$⁠. For a scalar-valued function $f(x,y,z)$ of multiple variables, $∇f$ refers to the gradient, while for a vector field $u→(x,y,z)$⁠, $∇•u→$ refers to the divergence. Further, let $N$ be a stoichiometric matrix, $v$ a vector of reaction rates and $s$ a vector of species concentrations. Equation (2) can be written in matrix form as

where the general rate vector, $v=v(s,p,r,t)$ is a function of species concentrations, a parameter vector (⁠$p$⁠), a position vector (⁠$r$⁠) and time. The vector $v$ can be partitioned into reaction equations, $ρ$⁠, and cross-membrane transport equations, $τ$⁠, while $s$ can be partitioned over $χ$ into sets of species grouped by compartment. Furthermore, $U=U(s)$ is a matrix with compartment-specific phase average velocities on the diagonal. Similarly, $D=D(s)$ is a matrix with diffusion coefficients on the diagonal. Each species has a unique diffusion coefficient.

Note that Equation (3) is a specific case of a more general conservation equation, which for an arbitrary variable, $φ$⁠, can be written as follows:

where $t$ is time and $∇$ is the gradient operator, $f$ and $g$ are functions of time, position (⁠$r$⁠), $φ$ and the gradient of $φ$ (⁠$∇φ$⁠); $f$ and $g$ are, respectively, called the flux and source.

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.

2.1.4 ‘Solute-driven advection’.

Long-distance photoassimilate transport in plants occurs via the phloem (Bieleski 2000, van Bel 2003). 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 (Bear 1972),

where $κ$ is the hydraulic conductivity (a measurable quantity assumed to be constant, Thompson and Holbrook 2003b), $μ$ 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 (Atkins and de Paula 2002),

where $Π$ is the osmotic pressure and $P0$ 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,

Furthermore, the contribution of the diffusion term to phloem translocation is typically two orders of magnitude smaller than the contribution of the advective term (Henton et al. 2002; Thompson and Holbrook 2003b), in other words

This allows phloem translocation to be modelled with the following approximation:

Thompson (2006) called this ‘solute-driven advection’.

The next section describes the equations governing catalytic conversions and cross-membrane transport, in other words, finding the rates to populate $v$ in terms of metabolite concentrations and parameters.

2.1.5 Reactions and cross-membrane transport.

To model biochemical reactions, we use a generic enzyme-kinetic rate equation, i.e. the generalized reversible Hill equation, which for a bi-substrate bi-product reaction reads (Hanekom et al. 2006; Rohwer et al. 2006, 2007)

where $s$⁠, $p$ and $m$ are, respectively, substrate, product and modifier concentrations scaled by their specific half-saturation constants. $Vf$ is the limiting forward reaction rate, $Γ$ is the mass-action ratio and $Keq$ is the equilibrium constant for the reaction (the disequilibrium ratio $Γ/Keq$ determines the direction of the reaction); $h$ is the Hill coefficient (Hofmeyr and Cornish-Bowden 1997; Hanekom 2006); $α$ and $γ$ are the scaling coefficients, either activating or inhibiting the enzyme through, respectively, a K-effect or V-effect (Cornish-Bowden 2012). The dot product, $m⋅r$⁠, is used for convenience.

If $α=γ=h=1$ and $m=0$ (i.e. there are no allosteric effects and the enzyme does not show cooperativity), then Equation (12) simplifies to

The equation for uni-uni reactions is still simpler:

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. $∂t∇e=0$⁠.

2.1.6 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

and appear as such in the model. Solute flow and reactions are defined on certain compartments and sometimes also span a number of compartments.

2.1.7 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 (Beek et al. 1999). 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 Henton et al. (2002). Similarly, Equation (11) is used as a simplified version of the conservation equation for modelling phloem (Thompson 2006).

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 generalized 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 (Hofmeyr and Cornish-Bowden 1997; Rohwer et al. 2006).

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 Bear (1972) to model phloem flow, and the generalized reversible Hill equation to calculate reaction rates.

2.2 Computational methods

Computational analyses were performed on a laptop computer using the Python programming language with the open-source numpy (Oliphant 2006), scipy (Virtanen et al. 2020) and FiPy (Guyer et al. 2009, 2017) libraries. FiPy is a PDE solver written in Python using finite volume methods, and provides an object-oriented interface to solving coupled sets of PDEs. Suppose a problem can be cast in the form of a general conservation equation,

where $ρ$ is some coefficient, $φ$ is the variable being solved, $u→$ is a vector coefficient (e.g. velocity), $Di$ is a generic coefficient of diffusion (e.g. 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 visualized with the matplotlib plotting library (Hunter 2007).

3. RESULTS

3.1 Example of building a model

Previous models of phloem flow have usually approximated the sieve tube as a cylinder with loading of solute at the one end and unloading at the other (Henton et al. 2002; Thompson and Holbrook 2003a). A few models have also considered branched vascular bundles (Lacointe and Minchin 2008). Some others have included general sink type demand, such as ‘respiration’ or ‘starch’ (Quereix et al. 2001). 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.

3.1.1 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 Fig. 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 remobilization. A futile cycle known to exist in sugarcane (Komor 1994) 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 synthesized 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 (M_sk/N_sk, representative of ATP/ADP; and H_sk/H_vc, referring to the movement of protons across the tonoplast, i.e. the membrane surrounding the vacuole). Importantly, the species H_sk is not meant to refer to protons in a mechanistic way, but is included to generically couple the uptake of S_sk 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 Fig. 2). A feedback mechanism is also present, reminiscent of the allosteric inhibition of phosphofructokinase by ATP (Evans et al. 1981). The import of S_sk into the vacule is via a proton antiporter (Reinders et al. 2006). The sink component of the reaction network in the model has a number of distinguishing features. The reactions occur in more than one compartment (Welbaum and Meinzer 1990; Beevers 1991; Rae et al. 2005a; Kruger et al. 2007), in contrast to previous models where only subcellular compartment was considered (Rohwer and Botha 2001; Uys et al. 2007). It is possible for some metabolite intermediates to enter futile cycles (Komor 1994; Zhu et al. 1997), while other metabolites may form part of moiety-conserved cycles (e.g. ATP/ADP) (Dancer et al. 1990; Cloutier et al. 2009). Assimilate can be partitioned between respiration or storage (Whittaker and Botha 1997; Bindon and Botha 2002), where accumulation is due to energy-driven antiporter proton gradients (Thom and Komor 1985; Reinders et al. 2006). Stored assimilate can also be remobilized from the vacuole (Moore 1995). 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 phosphofructokinase [PFK] and pyruvate kinase [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 side effect is to introduce an additional ATP-consuming or -producing reaction. This was the reason for including reaction 12.

3.1.2 Partial differential equations.

Table 2 summarizes 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.

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:

With the following shorthand,

the PDEs governing species behaviour in the symplast are:

Note that the stoichiometry for reaction 10 is $1P+2M→1Y+2N$ and for reaction 11 it is $1Y+4N→1Z+4M$⁠. Reaction 10 is representative of upper glycolysis, where two ATP are consumed, and reaction 11 of lower glycolysis where four ATP are produced. All other reactions have unitary stoichiometry. The specific case of assimilate translocation in the phloem is modelled by

The stoichiometric model in Fig. 2 has two conservation relations describing the total amounts of $Msk+Nsk$ and $Hvc+Hsk$⁠. These conservation relations were maintained with the following algebraic expressions:

The $vi$ appearing on the right-hand side of the PDEs are themselves functions and are determined by the rate equations (Equation (12)).

3.2 Discretization and parameter selection

3.2.1 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, $L=1.0$⁠, and divided into 50 finite volumes. The compartments, i.e. source (leaf), phloem, symplast and vacuole, were then defined on these finite volumes (Fig. 3). Each finite volume element has at least three 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 Fig. 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 10 elements to a node/internode pair; the first element of each group of 10 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 $N=50$ 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 nine times longer than a node, as approximation of real internode lengths (van Dillewijn 1952).

Because this was a core model, any units assigned to parameters and variables would be arbitrary to some extent. However, it should be emphasized 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 discretization 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.

3.2.2 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 Fig. 2, starting at the futile cycle in the model. Reaction 5, $S→P$⁠, is the opposite of the half reaction, $P→S$⁠, of enzyme 4. By definition, $KSP=KPS−1$⁠; therefore, $KSPKPS=1$⁠. It is possible for reaction 4 to have an equilibrium constant larger than $KPS$ because it is coupled to the half reaction, $M→N$⁠, which does have a large equilibrium constant so that $K4=KPSKMN$ can be large. Note, however, that the interconversion of M and N (in both directions) is involved in other reactions as well, which places a constraint on the choice of equilibrium constants for those other reactions as the equilibrium constant for this specific half reaction may not change.

The easiest way to satisfy all these criteria was to specify the equilibrium constants for the half reactions as follows: $KSP=KPS=1,KPY=0.1,KYZ=104$⁠, and $KMN=10$ (by implication $KNM=0.1$⁠). The high value of$KYZ$ allows reaction 11 to provide the driving force for the regeneration of the moiety $N→M$⁠. The equilibrium constant for the half reaction of transport of a single species across compartmental boundaries (i.e. not considering any energy coupling) is 1 by definition. With these choices, equilibrium constants for all reactions are equal to 1, except for reactions 4, 7, 10 and 12, which have a value of 10 (they are ‘energy-driven’ reactions).

3.2.3 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 five points along the stalk that produce assimilate through phloem loading, but 50 that consume it through phloem unloading. The chosen value was $Vf1,2=10$⁠.

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 remobilization 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 Fig. 4. Recall that $z$ is defined on the interval $z∈[0,L]$ and that $L=1$ in this geometry. The activities that increased with increasing $z$ are the phloem unloading steps and the vacuolar accumulation steps.

The rough guiding principle in selecting parameters was to observe accumulation of S_vc 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.

3.2.4 Flow coefficients.

Diffusion coefficients were set at the same value of $10−5$ for all species in the symplast. The flow coefficient was assigned a value of

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.

3.2.5 Boundary and initial conditions.

Boundary conditions were applied at $z={0,L}$ such that there is a zero gradient at the boundaries (e.g. for state variable $S$⁠, $∂zS(0)=0$ and $∂zS(L)=0$⁠), ensuring zero flux over the boundaries (the default boundary conditions in the FiPy software). Such boundary constraints are a special case of Neumann boundary conditions; see the Discussion for a distinction between Dirichlet (fixed value) and Neumann (fixed gradient) boundary conditions.

The initial values of each species variable, given in amounts, are summarized in Table 3. The amounts of the moiety-conserved species add up to a constant sum as follows:

3.3 Simulation results

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

Sawtooth profiles: The distributions of S_ph, S_sk, S_vc, P_sk and P_vc 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_vc, as can be seen on the left-hand side of Fig. 5. S_so concentrations were initially evenly distributed along the stalk, and end the simulation at a shallow gradient monotonically decreasing with $z$⁠. Similarly S_ph in general decreased down the stalk, except for the local peaks around the nodes. The S_sk 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_vcincreased along the stalk, with concentrations at immature end climbing faster than those at the mature end.

Moieties: The simulation profiles of moiety pairs, M_sk and N_sk as well as H_sk and H_vc, form mirror images across $z$ and $t$⁠. This is due to the moiety conservation constraints, i.e. the pairs always have to sum up to a constant (Equation 21).

Reaction rates: As with the concentrations, most of the reaction rates settled into stable profiles (at approximately $t>150$⁠). The rate of synthesis (⁠$v1$⁠) in the leaf increased with $z$ and $t$⁠; this was due to decreasing S_ph. Similarly phloem loading (⁠$v2$⁠) increased because of increasing S_so and decreasing S_ph. Phloem unloading (⁠$v3$⁠) fell to a level much lower than the loading rate, because only five source elements have to supply a total of 50 sink elements (recall that an element with a node defined on it overlaps a part of the symplast as there is parenchymal tissue at the nodes as well as in the internodal regions, see Fig. 3).

The rates of the cooperative, allosteric reactions $v10$ and $v11$ stabilized relatively quickly, reflecting the allosteric inhibition by M_sk—homoeostasis was almost immediate. Since $v11$ is responsible for regenerating the energy equivalents, homoeostasis will impact the vacuolar uptake step. Furthermore, the moiety regeneration and consuming reactions were very close to a steady state, i.e. $4v11−(2v4+v7+v10)≈_0$ for the entire length of the stalk at $t=600$ time steps.

The rate $v6$ 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 Fig. 8.

Futile cycling: The reactions $v4$ and $v5$ form a futile cycle in the symplast. The calculated rate of $v5$ was higher than $v4$⁠; therefore, more S_sk was being broken down than regenerated. Moreover, S is also broken down in the vacuole by $v6$⁠, adding to the overall degree of futile cycling as illustrated in Fig. 6. Note that in spite of this futile cycling ratio of $>1$ it was still possible for S_vc to increase over time, because S is additionally also taken up into the symplast through $v3$⁠.

3.4 The ‘sock experiment’

Increased photosynthetic rates in response to shading have been experimentally demonstrated in sugarcane (McCormick et al. 2006, 2008); shading of the leaves was argued to correspond to an increase in sink demand. All but one leaf of a sugarcane stalk was covered with shade cloth and the effect on various sugars and rates measured. Sucrose concentrations decreased in young internodal tissue, while the carboxylation efficiency and rate of electron transport increased in the exposed leaf (McCormick et al. 2006). Expression of a number of genes associated with photosynthesis, mitochondrial metabolism and sugar transport were also upregulated (McCormick et al. 2008).

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 $t=102$ during a time-dependent simulation of the reference model.

In other words, all but the second leaf had the maximal velocity for reaction 1 set to zero (Fig. 7). The following effects can be noted after the ‘shading’ event, compared to the reference model:

• $v1$ (sucrose synthesis in the leaf) and $v2$ (phloem loading) in the second (i.e. non-shaded) 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.

4. DISCUSSION

In this paper we have presented a general ADR 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 (Jensen et al. 2016), and recently has been extended to a whole-plant framework (termed CPlantBox, Zhou et al. 2020) 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.

4.1 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 (Thompson 2006). 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 (Atkins and de Paula 2002). A phenomenological equation has been proposed by Michel (1972) and is used by Thompson and Holbrook (2003a), Minchin and co-workers (Lacointe and Minchin 2008; Hall and Minchin 2013) and Seleznyova and Hanan (2018). Lacointe and Minchin (2008) 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.

4.2 Model

The model formulation in Fig. 3 is of course a simplification of what happens in a sugarcane plant, and the formulation could be criticized on this ground. In what follows, we provide a justification for our simplifying choices.

First, the intricate network of reactions in the leaf is abbreviated to one reaction in the model. This assumption is partly based on the finding that photosynthetic rates respond to changes in demand rather than the other way around (Wu and Birch 2007; McCormick et al. 2008).

Next, the model does not account for growth or volume changes with time. Large time scales therefore do not give an accurate picture of what happens in the plant. To model growth a method such as ‘adaptive meshing’ would be required to dynamically adapt the domain to the change in geometry (Plewa et al. 2005). Partial differential equations 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 $y=f(x,∂x)$ is some function of interest on the interval $[a,b]$⁠, then the Dirichlet boundary condition would be $y(a)=c1$ and $y(b)=c2$⁠, while the Neumann boundary condition is slightly more complicated, $∂xy(a)=c3$ and $∂xy(b)=c4$⁠. 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 discretization 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 (Geuzaine and Remacle 2009), a custom 3D finite element mesh generator, this comes of course at a vastly increased computational cost.

Finally, for a particular position $z$ on $L$⁠, $Hsk+Hvc$ is equal to the same constant as the same sum at an immediate, neighbouring position. Since H_sk is able to diffuse along the entire length of stalk, but H_vc is not (it is bounded by the vacuole), this means they do not necessarily have to sum to the same constant as one moves from element to element along the stalk. However, this is not actually a problem since the total amount of $H$ in the entire stalk does not change. Thus, there can in principle be a local violation of the moiety conservation (although we have not observed it in our data), but not across the length of the stalk.

To a lesser extent the same applies to M_sk and N_sk. 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.

4.3 Simulation

Figure 8 attempts to explain the cause of the observed parabolic reaction rate profiles. Suppose a reaction is sandwiched between two others, such as reaction 2 in the scheme on the left. If the maximal activity for the first reaction increases with length and that of the third reaction decreases, such as $V1$ and $V3$ in the plot on the right, then the highest rate that $v2$ can attain will be where $V1$ and $V3$ cross. Furthermore, if its own maximal activity (⁠$V2$⁠) remains constant along $z$⁠, close to the left-hand side the rate $v2$ cannot run faster than $v1$ without starting to deplete the linking metabolite/s, which will of course slow down $v2$⁠. By symmetry the same applies to $v2$ and $v3$⁠, and $v2$ is fastest somewhere in the middle of $L$ where it has a steady supply of substrate and its products are being removed fast enough not to become inhibiting. This could account for the observation that rates of growth and accumulation in medium mature tissue in sugarcane are the fastest, i.e. it could be a consequence of reciprocal profiles of up- and downregulation of enzymes along the stalk. An example of experimentally determined, steady-state fluxes that are higher in the medium mature internodes, rather than the surrounding internodes, is found in Bindon and Botha (2002): the flux to sucrose from the hexose phosphate pool was $2.85±0.40$⁠, $7.63±2.17$ and $4.73±1.12$ nmol hexose equivalents per mg protein per minute, for internodes 3, 6 and 9, representing immature, medium mature and mature tissue.

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 $v1$ in leaf 2, due to attenuation of $v1$ in the other leaves, is in keeping with the main experimental findings (McCormick et al. 2006, 2008), which showed an increase in photosynthetic rate (carbon dioxide assimilation and photosynthetic electron transport) in response to a sink perturbation. This can be directly explained by the quasi-linearity of the pathway. P_vc is accumulated at all sites along the domain, which leads to a local depletion of intermediates along the chain, specifically S_ph. This causes a decrease in the mass-action ratio of reaction 2 (⁠$Γ2$⁠), which increases the rate at which S_so is depleted and which then decreases $Γ1$⁠. The only place where this will have an effect is in leaf 2, as photosynthesis is completely blocked in the other leaves.

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_vc 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 Fig. 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.

5. CONCLUSIONS

We have presented a general ADR 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 (Uys et al. 2021).

SUPPORTING INFORMATION

Supplementary text document

• adr_supplementary.pdf Summary document with instructions for setting up the computational environment.

Code

• adr.py Model definition.

• run_model.py Script for model execution.

• LICENSE.txt Software licence.

Data

• adr_normal_data.npy Data for baseline model evaluation.

• adr_sock_data.npy Data for ‘sock experiment’ evaluation.

SOURCES OF FUNDING

This work is supported in part by funds from the South African National Research Foundation (NRF: # 81129 and # 114748) and the South African National Bioinformatics Network. The funding bodies played no role in the design of the study, collection, analysis and interpretation of data, decision to publish and in writing the manuscript.

ACKNOWLEDGEMENTS

The authors thank the anonymous reviewers for their valuable suggestions.

CONFLICT OF INTEREST

None declared.

CONTRIBUTIONS BY THE AUTHORS

L.U. and J.M.R. conceived the analysis, performed the simulations, analysed and visualized the results. J.M.R. procured funding. L.U. wrote the first draft of the manuscript. J.M.R. and J.-H.S.H. supervised the study. All authors reviewed and edited the manuscript.

DATA AVAILABILITY

Model code and simulation data are provided as Supporting Information.

This paper is dedicated to the memory of Prof Prieur du Plessis, whose help in conceiving and setting up this analysis has been invaluable.

LITERATURE CITED