Hydrodynamics of steady state phloem transport with radial leakage of solute

Long-distance phloem transport occurs under a pressure gradient generated by the osmotic exchange of water associated with solute exchange in source and sink regions. But these exchanges also occur along the pathway, and yet their physiological role has almost been ignored in mathematical models of phloem transport. Here we present a steady state model for transport phloem which allows solute leakage, based on the Navier-Stokes and convection-diffusion equations which describe fluid motion rigorously. Sieve tube membrane permeability Ps for passive solute exchange (and correspondingly, membrane reflection coefficient) influenced model results strongly, and had to lie in the bottom range of the values reported for plant cells for the results to be realistic. This smaller permeability reflects the efficient specialization of sieve tube elements, minimizing any diffusive solute loss favored by the large concentration difference across the sieve tube membrane. We also found there can be a specific reflection coefficient for which pressure profiles and sap velocities can both be similar to those predicted by the Hagen-Poiseuille equation for a completely impermeable tube.


INTRODUCTION
Phloem transport denotes long-distance transport, mainly of assimilates arising from photosynthesis, and is the movement of a solution in a continuum of interconnected cells, sieve elements, within the phloem of the vascular tissues in plants. It is currently accepted that solutes enter and exit the sieve tubes at sources and sinks, water enters and exits osmotically, and the solution moves in these sieve tubes due to the consequent osmotically generated pressure gradient: the theory of Münch pressure flow.
However, there is also a considerable radial exchange of solutes between the sieve tubes and the adjacent cells along the longdistance pathway between source and sink regions, the so-called transport phloem (i.e., main veins, petioles, stems, and main roots) (Van Bel, 2003). This radial exchange has hardly been addressed in mechanistic modeling of phloem transport. In early work when there were vigorous discussions on the pathway mechanisms the flux was specified (e.g., Tyree et al., 1974;Goeschl et al., 1976), and more recently models of cambial growth in trees have included radial solute flux as a variable (among many processes) (De Schepper and Steppe, 2010;Hölttä et al., 2010). The radial flux can take both apoplastic and symplastic pathways, and either path can predominate according to plant species and conditions (Patrick and Offler, 1996). Here we consider the apoplastic path only, where exchange is believed to be a leak/retrieval system (Tegeder et al., 2012). The passive leak, diffusive but possibly facilitated, is driven by the high concentration difference between the sieve element/companion cell complex and its surrounding apoplast Minchin and Thorpe, 1987;Patrick, 1990;Van Bel, 1990;Carpaneto et al., 2005;Thorpe et al., 2005). The active retrieval is driven by sucrose/proton symport (Hafke et al., 2013). There is good direct evidence for the mechanism of retrieval in several species but, for unloading to be studied, no technique has yet given appropriate access. Passive unloading of photosynthates into the apoplast has been estimated at about 6% cm −1 in bean (Phaseolus vulgaris) (Minchin and Thorpe, 1987;Van Bel, 1990). This passive radial exchange of solutes is coupled with radial water flux. For example, leakage of radio-labeled photosynthates responds to changes in the apoplastic water potential [e.g., by perfusing the apoplast with sorbitol ; mannitol (Aloni et al., 1986;Cabrita, 2011) or polyethylene glycol (Cabrita, 2011)]. The leakage of other metabolites shows similar behavior (Aloni et al., 1986).
Given the current information about the complex process of solute exchange in the transport phloem, and in the belief that any realistic model for transport phloem should include such a process, we chose as a first step to consider only the passive leak, characterizing the sieve tube membrane by a membrane permeability P s (in our case for sugars) and a reflection coefficient σ which describes the extent to which solute exchange is affected by the concentration difference across the membrane (Nobel, 2009). In all existing mathematical models of long-distance phloem transport based on the Münch pressure flow hypothesis [see Table  1 of Thompson and Holbrook (2003a)] the reflection coefficient was set to σ = 1, thus completely decoupling water exchange and solute exchange by assuming an ideally semipermeable membrane. The omission of radial solute exchange is not only a major gap in the description of phloem transport, but also makes the existing mathematical models suspect for interpreting the vast amount of experimental data that has been collected concerning phloem physiology.
Here we investigate the role of radial solute exchange in the transport phloem by studying a relatively simple system: steady state transport of a homogeneous solution (consisting of a single solute and water) in a tube surrounded by a bath of solution with predefined (but not necessarily constant) pressure and concentration. The simplicity of this setup allows us to highlight the effect of including radial solute exchange by comparing our findings with previous results from Thompson and Holbrook (2003b) and Phillips and Dungan (1993). In order to model the hydrodynamics of this system in a rigorous way we adapted a formalism introduced by Phillips and Dungan (1993), based on the Navier-Stokes and convection-diffusion equations. We extended this approach by introducing boundary conditions that specify radial exchange of solute across the sieve tube membrane. These exchange processes were described using the formalism of irreversible thermodynamics for transport across membranes (Kedem and Katchalsky, 1958). Our main focus is to study the importance for phloem transport of water and solute movement through the sieve tube membrane, and to evaluate the common practice in phloem flow modeling to assume semipermeable tubes.

MODEL ASSUMPTIONS
We propose a hydrodynamic model of phloem transport, according to the Münch pressure flow hypothesis, with allowance for water and solute exchange along the pathway corresponding to the transport phloem (Van Bel, 2003) (Figure 1). The following assumptions are made: 1. The sieve tube is considered as a right circular cylinder of length L and radius R, such that R L, limited by a porous membrane through which water and solutes fluxes occur. We use cylindrical coordinates where the axis along the pathway is denoted by z = [0, L], the radial coordinate is r = [0, R] and ϕ = [0, 2π] is the azimuth. The tube is surrounded FIGURE 1 | Sieve tube model. Solution enters the sieve tube with velocity U, pressure p i and concentration c i . There is radial exchange of water and solute across the sieve tube membrane along the pathway. The sieve tube is surrounded by a medium ("apoplast") with pressure p out (z) and concentration c out (z). R and L are sieve tube radius and length, respectively. r and z are the directions of radial and axial flow, respectively. by unspecified tissue (for simplicity called apoplast in the following). 2. The flow is axisymmetrical, i.e., flow depends on axis coordinate z and radial coordinate r, but not on azimuth ϕ. 3. Sieve tube sap is regarded as a homogeneous solution of a single solute with concentration c(z, r) in water and assumed to be an incompressible Newtonian fluid of constant density, ρ.
Sieve tube sap viscosity, μ, is constant, i.e., independent of solute concentration. 4. End effects caused by the entry and exit of sap in the sieve tube are negligible. 5. The system is at steady state, i.e., time independent. 6. The sap enters the sieve tube with an average speed U, with an average solute concentration c i and at average turgor pressure p i . 7. The membrane hydraulic conductivity, L p , is constant. Solute exchange across the membrane is regarded as a passive process in which there is: (i) diffusion of solutes that is linear with the concentration difference between the sieve tube and the surrounding apoplast, and dependent on the sieve tube membrane solute permeability, P s ; (ii) convection of solutes dragged by water through the membrane. 8. Apoplast solute concentration, c out (z), and hydrostatic pressure p out (z) are not affected by solute or water exchange across the sieve tube membrane. These functions have to be provided as boundary conditions. 9. There is a constant water potential gradient d out dz surrounding the sieve tube. 10. Diffusion of solute within sieve tube sap is isotropic and obeys Fick's law of diffusion with a constant diffusion coefficient D, solute-specific. 11. There is no slip at the sieve tube membrane. 12. There are no chemical reactions. 13. For simplicity, we assume osmotic pressure to be given by van't Hoff equation for dilute solutions: (z) = R g Tc (z), where R g is the universal gas constant and T the absolute temperature. 14. Sieve plates are transverse and spaced at regular intervals (Figure 1). The effect of sieve plates on the overall sieve tube conductance is described by an impedance factor β which is related to sieve tube element structure and given by Thompson and Holbrook (2003a): where α is the fraction of sieve plate area composed of sieve pore area; l is the length of sieve tube elements in between sieve plates; r p and l p define the sieve plate pore dimensions: radius and length respectively (Thompson and Holbrook, 2003a). As α = 1, l > l p and R > r p we have always β < 1. For the special case of no sieve plates the impedance factor is set to β = 1.

FUNDAMENTAL EQUATIONS
As in Phillips and Dungan (1993) the calculation starts from the steady state Navier-Stokes equation for an incompressible Newtonian fluid which, using the assumption of axisymmetrical flow, can be written in cylindrical coordinates as (Truskey et al., 2009) ρ u r ∂u r ∂r where pressure p and velocities u r and u z are functions of axial coordinate z and radial coordinate r. In most cases we omit these dependences in our notation for simplicity. The effect of sieve plates is described introducing an effective viscosity μ * = μ/β, thus having a more viscous fluid as β = 1 [see Equation (1)]. This approach treats the sieve tube as a conduit of essentially uniform resistivity, a close approximation to reality in view of the distribution of the sieve plates within sieve tubes (Weir, 1981;Henton et al., 2002). Pressure p includes the gravitational effect, i.e., p = p + ρgh; where p' is termed the hydrodynamic pressure inside the sieve tube, g the local acceleration of gravity and h is the vertical coordinate above a standard reference plane. Note that in contrast to Phillips and Dungan (1993) we use the Navier-Stokes equation including inertial forces [left hand side of Equations (2) and (3), respectively] to keep the formulation as general as possible. The actual impact of inertial forces is discussed in the Results section. The continuity equation for incompressible fluids in cylindrical coordinates is (Truskey et al., 2009) 1 r and the convection-diffusion equation for dilute solutions in steady state is (Truskey et al., 2009) where concentration c is a function of coordinates z and r in the same way as pressure and velocities. The axial and radial fluxes of solute inside the sieve tube are defined as sum of convective and diffusive components

BOUNDARY CONDITIONS
To allow for solute and water exchanges across the sieve tube membrane, the following boundary conditions are used: radial flux of solution through the membrane (at r = R) is given by Starling's equation (Nobel, 2009): where the van't Hoff expression of osmotic pressure for dilute solution has been inserted. The reflection coefficient of the sieve tube membrane σ assumes a value between 0 (totally permeable) and 1 (semipermeable). Pressure, p out , and solute concentration c out , outside the tube do not to depend on r and are assumed to vary linearly with distance: and Radial passive flux of solutes through the membrane on the one hand is given by Equation (7) for r = R, on the other hand it can be defined as the sum of convective flux (solute dragged by solvent) which is linear with the average solute concentration in the membrane,c(z), plus diffusive flux across the membrane that is described by the sieve tube membrane solute permeability, P s , and the concentration difference across the membrane (Kedem and Katchalsky, 1958;Benedek and Villars, 2000).
The concentration of solute,c(z), within the sieve tube membrane is given byc which means that the fluid velocity will be in the radial direction only, at the boundary. Symmetry at the center of the tube implies that radial velocity as well as radial derivatives of axial velocity and concentration are zero at r = 0: u r (0, z) = 0, ∂u z ∂r (0, z) = 0 and ∂c ∂r (0, z) = 0. According to our model assumptions, flow is already fully developed at the beginning of the transport phloem region considered in our model; loading of solutes and water has occurred in the source region, not being part of our model upstream. Therefore, the boundary conditions at the inlet of the tube z = 0 are u z (0) = U, c(0) = c i and p(0) = p i where the bracket denotes the average over the cross section c(z) = 2 R 2 R 0 c(r, z)rdr. Note that the hydrodynamic equations (2) to (4) describing the solution movement inside the sieve tube do not depend on the solute concentration explicitly, and we take viscosity to be independent of concentration. The flow depends on concentration through the boundary condition Equation (8) only. Thus, we can treat independently the hydrodynamics problem that describes how pressure and velocity change, and the solute transport equation (5) which is solved once the velocity field is determined.

DIMENSIONAL ANALYSIS
As with many problems in fluid dynamics and especially when applying the Navier-Stokes equation, the problem becomes more simple to solve using dimensional analysis (Regirer, 1960;Kundu and Cohen, 2008). With this method, it is possible to predict physical parameters that influence the sap flow and determine the relationships between several variables (pressure, velocity and concentration) when an exact functional relationship is unknown. This is not possible with a direct numerical solution of the governing equations (2) to (5). Following Regirer (1960), Phillips and Dungan (1993) and Thompson and Holbrook (2003b), we use system geometry and boundary conditions as scales to define dimensionless variables in the following way: Inserting these definitions into the governing equations, the corresponding scaling of the remaining variables and parameters [unique except for numerical prefactors which were chosen according to Thompson and Holbrook (2003a)] becomeŝ Substituting the dimensionless variables into the governing equations, the dimensionless Navier-Stokes Equations (2) and (3) become: Re ε û r ∂û z ∂r +û z ∂û z ∂ẑ Here ε = R/L is an abbreviation for the aspect ratio and Re = ρRU/μ * is the Reynolds number, the ratio of inertial to viscous forces. The continuity equation (4) now reads 1 r The convection diffusion equation (5) becomes where Pe r = RU/D is the Péclet number for radial flow. The components of the sieve tube solute flux are given by: where Pe z = LU/D = Pe r /ε is the Péclet number for axial flow. The Péclet number for any direction gives the ratio of the rate of convection of the solute in that direction by the flow of sap, to the corresponding rate of diffusion of that solute driven by a concentration gradient.
In dimensionless notation, the boundary conditions Equations (8), (11) and (12) at the sieve tube membrane arê The dimensionless parameterĤ = 8L Uμ * , defined byˆ =Ĥĉ, gives the ratio of osmotic forces to viscous forces for a flow with a characteristic velocity U (Phillips and Dungan, 1993).

SERIES EXPANSION
Proceeding in the same way as Phillips and Dungan (1993) we expand the state variables,û r ,û z ,ĉ andp as power series of the small dimensionless aspect ratio ε: The accuracy of the expansions, i.e., the number of terms to include, will depend on the value of ε: the smaller ε, the more significant are the first terms compared to higher order terms. In the case of phloem transport, ε is smaller than 10 −3 for typical sieve tube dimensions (Esau, 1969). For this reason we will only consider the first two terms (zeroth and first order) of the expansion to describe flow in a sieve tube. Higher order terms might be of interest when studying effects on a smaller length scale (e.g., single sieve tube elements).
Collecting terms of the appropriate order (ε 0 and up to ε 1 , respectively), using respective boundary conditions and inserting zeroth order results into the calculation of first order leads to the following expressions: axial velocity can be calculated from Equations (15) and (22) aŝ (24) where pressuresp 0 andp 1 are independent of radial coordinate r as a result of Equation (14). Using Equations (23) and (24), radial velocity follows from Equation (16) aŝ Inserting Equation (20) into Equations (25) and (26) the relation of pressure and concentration becomes The average of zeroth order axial velocity over the cross section is û z0 = − ∂p 0 ∂ẑ . Thus in our model the pressure gradient at any distance from the origin is linear with the average axial velocity, but dependent both explicitly and implicitly on the radial flux of water at that location, given by the zeroth order functionû lat 0 ẑ , see Equation (27).
Evaluating the convection diffusion equation (17) for ε = 0 together with the boundary condition, we concludeĉ 0 ẑ to be independent of radial coordinate: ∂ĉ 0 ∂r = 0. Using this result and Equation (23) the first order expression of Equation (17) becomes The first order expression of the boundary condition Equation (21) is Combining Equations (29) and (30) and inserting Equations (25) and (27) leads to the following zeroth order equation for concentration The final first order expression for concentration follows from Equation (29) using the boundary conditions ĉ 1 (z = 0) = 0 and As the velocity profile depends on pressure only [Equations (23) to (26)] and first order concentration follows from zeroth order expressions [Equation (32)] we need to knowĉ 0 ,p 0 , andp 1 only in order to obtain the pressure and velocity profiles and thus describe phloem flow up to first order of the expansion in aspect ratio ε. These variables are obtained by solving the system of coupled differential equations (27), (28) and (31). For a semipermeable membrane, i.e., σ = 1 andP s = 0, the zeroth order equations (23), (25), (27) and (31) are identical to findings of Phillips and Dungan (1993) except for different numerical prefactors stemming from slightly different definitions of the dimensionless quantities. Also in this case Equations (27) and (31) are equivalent to the governing equations of Thompson and Holbrook (2003b), the only difference being in the boundary conditions. In this light it is not surprising that the zeroth order results Equations (27) and (31) for a permeable membrane can also be obtained using an approach similar to Thompson and Holbrook (2003a), based on conservation equations and a local application of the Hagen-Poiseuille equation together with appropriate boundary conditions, see Appendix B.

NUMERICAL ANALYSIS
As there is no analytical solution for the system of differential Equations (27), (28) and (29), they were transformed into a first order system of differential equations [Appendix C, Equations (C.3)-(C.7)] and solved with MATLAB (R2010b, MathWorks, Inc.) ode15s differential equations solver routine. The calculated zeroth order of concentration,ĉ 0 , turgor pressure,p 0 , the first order of turgor pressure,p 1 , and their respective derivatives were then used to determine the average first order concentration ĉ 1 , Equation (C.8). Concentration profile and sieve tube turgor pressure were determined as sum of zeroth and first order. From these profiles, average over cross section of axial velocity Equation (C.9), radial flow at sieve tube membrane Equation (C.10), axial solute flux Equation (19) and radial solute flux Equation (18) were determined.

PARAMETER VALUES
Unless otherwise specified, the parameter values of Table 1 were used (the basis for these values is given in Appendix A). The values of the non-dimensional model parameters related to the parameter values of Table 1 are shown in Table 2. Figure 2 shows profiles of pressure, concentration and flow for different values of the reflection coefficient σ. In the case of a semipermeable membrane (σ = 1), concentration and pressure differences across the sieve tube membrane (c and p are higher than c out and p out ) create a water potential difference, , across the sieve tube membrane that always draws water into the sieve tube (û r < 0, Figure 2D), causing a pronounced non-linearity of pressure in the direction of flow (Figure 2A). In the case of a membrane that is permeable to solutes as well as water (σ < 1), the water influx depends on the value of the reflection coefficient σ and can even be zero (Figure 2D), leading to a pressure profile almost identical to Poiseuille flow (Figure 2A). Boundary condition Equation (20) shows that a water influx (û r < 0) will occur as long as:

PERMEABLE MEMBRANE-THE EFFECT OF RADIAL SOLUTE EXCHANGE
For typical physiological conditions (Table 1), i.e., for the order of magnitude of pressure and concentration expected for sieve tubes, water influx into the sieve tube occurs only if σ > 0.7 ( Figure 2D). If the permeability of the membrane, P s , increases www.frontiersin.org December 2013 | Volume 4 | Article 531 | 5 Universal gas constant, R g J K −1 mol −1 8.314 Water density, ρ w kg m −3 998 a Thompson and Holbrook (2003a). b Phillips and Dungan (1993).
c Eszterle (1993). (with σ decreasing from σ = 1 to σ = 0.7; see Appendix A for the relationship of P s and σ), the absolute value of the pressure gradient in the direction of flow becomes smaller with decreasing axial velocity ( Figure 2C) as û z ≈ − ∂p ∂ẑ , see Equation (C.9). This trend arises from two factors that decrease the sieve tube solute concentration along the axis: first, the dilution created by water influx, dependent on the water potential difference across the sieve tube membrane; second, the passive efflux of solutes across the sieve tube membrane, dependent on the concentration difference between the sieve tube and the apoplast. The passive loss of solutes is favored by a higher sieve tube solute concentration compared with the surrounding apoplast. The decrease in concentration, as one moves further down the tube, means less water will enter because the water potential difference across the sieve tube membrane decreases, and more so for a more permeable membrane ( Figure 2D). Consequently, due to the volume conservation, the increase in axial velocity is less for more permeable membranes ( Figure 2C) and there is also a decrease in axial solute flux ( Figure 2E). Due to the concentration difference across the sieve tube membrane there is solute efflux which increases for smaller values of the reflection coefficient, σ, corresponding to a leakier membrane ( Figure 2F). However, the sieve tube membrane solute permeability, P s , is in the order of 10 −10 m s −1 (see Appendix A), and too small to cause dramatic changes in the sieve tube concentration. If the permeability of the sieve tube membrane, P s , increases (with a corresponding decrease in the reflection coefficient, σ) there is a smaller osmotic effect on water exchange, leading to less water entering the sieve tube ( Figure 2D) and a smaller decrease of solute concentration ( Figure 2B). Eventually, for a very permeable membrane there is a reversal and water moves out radially for σ < 0.7 (Figure 2D), and the axial velocity decreases along the axis (Figure 2C). Simultaneously, the pressure inside the sieve tube decreases less and less for smaller values of σ (Figure 2A). Axial velocity and axial solute flux at one point become null (Figures 2C,E). Numerical calculation has to stop at this point since, without specifying boundary conditions at the outlet of the tube, it cannot be decided whether there is inflow from the opposite side or the axial velocity remains zero along the tube.
The expected behavior of solute concentration ĉ for very permeable membranes (small values of σ) would be that ĉ decreases less with distance, approaching a value constant over distance for σ = 0, but not increasing over the initial value at z = 0. Obviously our numerical results of solute concentration ( Figure 2B) deviate from this expected behavior. The steep increase in ĉ ( Figure 2B) for very small values of σ stems from the fact that the pressure gradient on the left hand side of Equation (31) approaches zero much faster than radial solute flow on the right hand side. The extreme of this scenario is seen for a totally permeable membrane (σ = 0), in which case the radial water exchange is driven solely by the pressure difference across the sieve tube membrane, and the solute flux occurs predominantly through convection by the water moving across the sieve tube membrane. In this case Equation (31) becomes ∂p 0 ∂ẑ ∂ĉ 0 ∂ẑ = −û lat 0 /2 +P s ĉ 0 −ĉ out , concentration inside the sieve tube c 0 (z) should stay constant, and apoplast solute concentration c out (z) would be expected to tend to c 0 (z) as a consequence of the high permeability of the membrane. However, our model assumes c out (z) to be independent of radial solute flux, with constant apoplastic concentration gradient over distance [model assumptions 8 and 9, see also Equation (10)], leading to a higher concentration difference between inside and outside than expected, and, as a mathematical consequence of Equation (31), an accumulation of solute inside the tube with distance. Thus the results for solute concentration ( Figure 2B) and radial solute flow (Figure 2F), though consistent with the model assumptions, are not realistic for small values of σ.

COMPARISON OF ZEROTH AND FIRST ORDER OF EXPANSION
In order to justify the truncation of the series expansion we compare zeroth and first order results. For velocity and pressure Equations (23) to (28) the first two orders structurally only differ by expressions proportional to the Reynolds number Re which is quite small in our case (Re = 8.94 × 10 −5 for the values reported in Table 1). Thus we expect terms proportional to Re not to contribute significantly for the set of parameters given ( Table 1). The remaining terms of first order in Equations (24), (26) and (28) should be of the same order of magnitude as zeroth order, depending on boundary conditions. Together with the fact that the aspect ratio ε is small (ε = 2 × 10 −5 for the values of Table 1) we conclude that zeroth order strongly dominates the series expansions of velocity and pressure profiles. The comparison is less obvious for the concentration profile where first order Equation (32) is proportional to the Péclet number Pe r . Since the Péclet numbers in our case are Pe r = 3.62 and Pe z = 1.81 × 10 5 , it is obvious from Equation (19) that Table 1 were used for the numerical calculation. The diagonal indicates same order of magnitude of zeroth and first order coefficients in the power series expansion.
Clearly, diffusive processes (in this case passive) can be important for solute exchange between the sieve tube and its surroundings, and especially if the radial convection flow is small. Figure 3 shows a comparison of zeroth and first order concentration, i.e., the numerical solutions of Equations (31) and (C.8) for the parameter values of Table 1 and different values of the reflection coefficient σ. First order concentration never exceeds the order of magnitude of zeroth order. Together with the small aspect ratio ε it is thus clear that first order in general does not contribute significantly for this set of parameters (Table 1). This implies that out of the set of dimensionless parameters governing the model behavior together with the three boundary condition valuesp i ,p out ẑ andĉ out ẑ , there are only three independent significant ones:L p , H and the reflection coefficient σ. PermeabilityP s depends on σ, see Equation (A.2). Aspect ratio ε, Reynolds number Re (i.e., the effect of inertial forces) and Peclét number Pe r play a role only when the first order of the series expansion becomes relevant, such as in wider sieve tubes for which ε is bigger than we have considered (Table A1).

SENSITIVITY ANALYSIS
For a semipermeable membrane (σ = 1), the exploration of the model behavior for a wide range of possible values of the dimensionless parameters has been provided by Thompson and Holbrook (2003b). Extending this approach to semipermeable membranes would be laborious due to the additional parameters involved. Instead we performed a sensitivity analysis of our model in order to compare the significance of the model parameters, with particular focus on the role the reflection coefficient plays compared to the other parameters, varying each model parameter and boundary condition stepwise from −50 to 50% in relation to the parameter values of Table 1 [which are somewhat similar to the standard parameter set of Thompson and Holbrook (2003b)]. Figure 4 shows the relative change of the state variables pressure, average solute concentration, average axial velocity and radial velocity at the sieve tube membrane. Variation of the parameters causes the dimensionless values of pressure at the inlet of the tubep(z = 0) to differ even if initial pressure p i is kept constant. To avoid this effect we converted dimensionless results back into original dimensions for the sensitivity analysis. Since the values of the state variables (except for radial velocity) are fixed by boundary conditions at the inlet of the tube, see Figure 2, changes at the end of the tube (i.e., z = L) are taken as indicator of model sensitivity. Boundary conditions p i , c i , and U correspond to state variables p, c , and u z , respectively. In cases where these boundary conditions vary, curves of related state variables are not shown in Figure 4 because the variation of initial values dominates the changes of these variables. Variation of parameters and boundary conditions is straightforward except for the reflection coefficient σ ∈ [0, 1] which is the only parameter limited by a lower as well as an upper boundary. We chose σ = 1 as reference and show negative variation only (such that −50% corresponds to σ = 0.5). Doing so, σ overall appears as one of the most sensitive model parameter, see Figure 4. Given the difficulty of comparing results for varying σ with other (only partially limited) parameters, this result should be treated cautiously regarding a direct numerical comparison of the fractional sensitivities. Still it shows that even small variations of sieve tube solute permeability have a strong effect on the resulting system behavior. Thus there remains no reason to ignore radial solute exchange as it apparently plays in the same league as hydraulic conductivity, geometry and other system properties.
Another obvious result from Figure 4 is that apoplast pressure, concentration and respective gradients ( Table 1) do not contribute to the system behavior as strongly as other parameters and boundary conditions. The apoplast surrounding the phloem has been considered to have a constant water potential in most phloem transport models; very few authors (Tyree et al., 1974;Weir, 1981;Thompson and Holbrook, 2003b) considered the more realistic situation of a water potential varying with height. Our result suggests that a water potential gradient can be disregarded, since it constitutes the least influential parameter for all state variables. Nevertheless, in a branched architecture the apoplastic water potential of alternative sinks may well be important (Kaufman and Kramer, 1967;Lang and Düring, 1991).

DISCUSSION
Solute exchange across sieve tube membranes in the pathway region has hardly been considered in mathematical modeling of phloem transport. Tyree et al. (1974) and more recently Hölttä et al. (2010) and De Schepper and Steppe (2010) have considered a radial solute flow, but in each case the reflection coefficient's value was unity, thus both over-estimating the osmotic contribution to the driving force for solution flow (Equation 8), and decoupling solute and water fluxes (Equation 11). In this study we present a rigorous way of including radial fluxes of both water and solute in the pathway region, in order to give a more realistic analysis of phloem transport. As a first step, we consider passive solute exchange, as a sum of diffusive and convective fluxes across the sieve tube membrane, see Equation (11). Our model considers a diffusive solute flux, described by the irreversible thermodynamics formalism which relates membrane solute permeability, P s , hydraulic conductivity, L p , and solute reflection coefficient, σ (Appendix A) (Kedem and Katchalsky, 1958). In the following we compare our approach to the two main papers it is an extension of, i.e., Phillips and Dungan (1993) and Thompson and Holbrook (2003b). In both papers there are just two dimensionless parameters governing the system behavior,L p andĤ [R andF in the notation of Thompson and Holbrook (2003b)], whereas we have three independent parametersL p , H, σ when calculating the zeroth order of the series expansion, together with three boundary condition valuesp i ,p out ẑ and c out ẑ . Aspect ratio ε, Reynolds number Re and Peclét number Pe r enter the final equations if calculating the first order of the series expansion, as we did in contrast to Phillips and Dungan (1993). The main difference to Thompson and Holbrook (2003b) (except for allowance of radial solute exchange) is in the boundary conditions: with a permeable membrane, solute flux density is no longer constant along the transport pathway, thus an assumption of a fixed concentration value at the end of the tube becomes numerically intricate when comparing different values of membrane permeability. We therefore chose the boundary conditions of Phillips and Dungan (1993), i.e., fixed values of flux, concentration and pressure at the inlet of the tube. Unlike Phillips and Dungan (1993) and Thompson and Holbrook (2003b) we refrained from exploring the full range of possible values of the dimensionless parameters but focused instead on the effect of different permeabilities with their different rates of radial solute exchange, given a typical and representative set of physiological parameters.
For a semipermeable membrane (σ = 1), radial water exchange across the sieve tube membrane leads to a strong deviation from Poiseuille flow (Figures 2A,C), the more so the longer the sieve tube. For very long sieve tubes the model predicts that pressure can become negative [ Figure 5A, this phenomenon was called "mathematical plasmolysis" by Thompson and Holbrook (2003b)] and an unrealistically steep increase of axial velocity ["runaway phenomenon", Tyree et al. (1974)] is observed ( Figure 5C). These effects are weakened for the more realistic case of a permeable sieve tube membrane (σ < 1), see Figure 2. According to Equation (33), for a certain value of the reflection coefficient, i.e., for σ ≈ 0.73 for our choice of parameter values (Table 1), there is almost no water influx [û r r = 1 ≈ 0, see Figure 5D], pressure stays positive ( Figure 5A) and the "runaway phenomenon" of axial velocity does not occur ( Figure 5C). These results show that longer sieve tubes can sustain more stable flow conditions than sieve tubes limited by semipermeable membranes if solute efflux occurs, because of the smaller radial water flux (Figures 2D and 5D) which leads to a less severe dilution of sap (Figures 2B and 5B). At the special value of solute permeability corresponding to σ ≈ 0.73, pressure and axial velocity get closest to Poiseuille flow, and there is (approximate) radial water potential equilibrium as can be seen from Equation (20).
The assumption that sieve tubes are in water potential equilibrium with their surrounding apoplast in plant stems is generally made when measuring turgor pressure gradients indirectly for sieve tubes (Kaufman and Kramer, 1967;Wright and Fisher, 1980;Sovonick-Dunford et al., 1981). Murphy (1989) showed that symplastic connectivity between sieve tubes and the surrounding tissue would not greatly affect the approximation of radial water potential equilibrium. Although our results refer to an idealized sieve tube with no active solute transport, we conclude that this common assumption is more realistic in phloem transport modeling if sieve tubes are taken as permeable, because the value of the water potential difference, , across the sieve tube membrane is approaching zero only for a permeable membrane (Figure 5D), so that ≈ 0 for σ ≈ 0.73. Assuming this quasi radial water potential equilibrium, and since =out , the sieve tube axial turgor pressure gradient is simply given by the sum of the axial gradients, of the external water potential and the internal osmotic pressure. That is: Ignoring between sieve tubes and the apoplast, Equation (34) should yield good estimates of phloem turgor pressure and turgor pressure gradients specifically in pathway regions, as in stems.
Clearly, metabolism and storage processes occur throughout the stems of plants, especially in big specimens like trees, and so significant radial solute fluxes must occur (Minchin and Thorpe, 1987;Van Bel, 1990;Hoch et al., 2003). Our value of sieve tube membrane permeability to solutes, P s , (in this case sugars) lies in the bottom range of the values reported for plant cells, 10 −10 to 10 −6 m s −1 , and it is smaller than the values reported for noncharged solutes (Diamond and Wright, 1969;Nobel, 2009). Our choice of parameter values is justified in Appendix A. This small sieve tube membrane permeability to solutes reflects the efficient specialization of sieve tube elements to bulk flow, as it minimizes the diffusive solute loss across the sieve tube membrane. Hence, we have a more efficient system of carbon transport along the plant body that is built to keep solute losses at minimum. Such an efficient transport system is especially important for big specimens where source and sink regions are far apart. For example, in the case of Cucurbitaceae, the combination of a small solute sieve tube membrane permeability to sugars, e.g., sucrose, with the fact that they transport oligosaccharides that exist in the phloem only and seem not to leave it diffusively Gorham, 1964, 1965;Schaffer et al., 1996), further illustrates a specialization of sieve tubes to a more efficient system for solute transport bulk flow.

Frontiers in Plant Science | Plant Biophysics and Modeling
December 2013 | Volume 4 | Article 531 | 10

OUTLOOK
In this study we have applied the Navier-Stokes equation with boundary conditions that describe radial fluxes of both water and solutes to allow a better understanding of phloem transport dynamics. Reloading of solutes through active processes could be added to the boundary condition Equation (11). Whether described by the product of a transport coefficient and the apoplast solute concentration, or by Michaelis-Menten kinetics, the effect of solute reloading on phloem flow would be mainly to counteract the passive solute loss and draw less water into the sieve tube ( Figure 2D). Consequently the changes of turgor pressure, velocity and concentration with axial distance would be less than our predictions where there is only solute loss. Also, it would further illustrate the benefit for bulk flow in sieve tubes from minimizing radial solute loss. Thus the changes seen in both the axial velocity and axial solute flux would be attenuated and their profiles would be closer to a more stable system (Figures 2C,F). Additionally, if the reloading of solutes were taken into account, it would be necessary to formulate mechanisms that affect apoplastic concentration in relation to sieve tube solute concentration. This also would help to obtain a more realistic model behavior for very permeable membranes (σ < 0.5). By testing several hypothetical reloading mechanisms (e.g., linear function, Michaelis-Menten) together with possible functions that describe apoplastic concentration over distance, we could gain insight into how plants may have specialized in order to cope with scenarios or situations in which their membrane permeability is affected. The resistance imposed by sub cellular structures, i.e., cytoplasmic components, phloem proteins (P-protein), and their relation to sieve plates, could be added to the sieve plate impedance factor β, affecting the sieve tube resistance only, and described by the viscous term of the Navier-Stokes equation. The effects observed in our model would thus be enhanced, i.e., the pressure changes due to viscous losses. Hence, the problem turns out to be how to best describe the influence of those sub cellular structures on the overall resistance of the pathway. It would also be informative to investigate the effect of a different sieve tube element structure e.g., where surface arrangement affects the areas for flows in membrane and lumen, as with inclined sieve plates.
For the interpretation of phloem flow measurements, it is of significance that our model shows that, out of all the parameters, a crucial one, the reflection coefficient, typically is the least well known. Its importance suggests that it might be possible to estimate the value of the reflection coefficient by calibrating the model with experimental data, if other relevant parameters are sufficiently well known.

APPENDIX A. PHYSIOLOGICAL PARAMETERS SIEVE TUBE STRUCTURE
Sieve elements are elongated with the longitudinal axis parallel to the bundle of vascular tissues. Their length is several to many times greater than their width, but both dimensions can vary considerably within the same plant and between species and genera. Typically, the sieve element length varies from 100 to 5000 μm while its width varies between 10 and 400 μm, for an approximately cylindrical shape (Esau, 1969;Parthasarathy, 1975). Therefore, for a sieve tube length, L, of 0.5 m, comparable to a small plant, we specified the following dimensions of the sieve tube elements: radius R = l0 μm; length l = 250 μm; with 0.5 as the fraction of sieve plate area composed of sieve plate pore area; sieve plate pore radius r p = 0.23 μm; sieve plate pore length, which equals sieve plate thickness, l p = 0.5 μm (Thompson and Holbrook, 2003a,b). With these dimensions the sieve plate impedance factor β [Equation (1)] is 0.079, which accounts for the sieve plate contribution to the total sieve tube axial conductivity.

SAP VISCOSITY
The phloem sap viscosity depends on both temperature and phloem sap chemical composition (proteins, sugars and other solutes). According to Chirife and Buera (1997) the viscosity of sugar solutions at a given temperature T is given by: where μ 0 (T) is the viscosity of the solvent, in this case water, at absolute temperature T; x s is the mole fraction of sugars and a(T) and E(T) are parameters depending on both temperature and sugar species. For a temperature of 22 • C (295.15 K), μ 0 = 9.548 × 10 −4 Pa s (Kundu and Cohen, 2008), a = 0.905 and E = 57.19 (Chirife and Buera, 1997). However, for convenience, as mentioned in the model assumptions we consider phloem sap viscosity as constant (μ = 1.5 mPa s) well inside the range (1-2 mPa s) for typical values of sucrose concentration in sieve tubes, 300-900 mol m −3 (Taiz and Zeiger, 1998).

TURGOR PRESSURE
There are not many published results on phloem turgor pressure. Considering both the nature and location of phloem in the plant body, it is clearly difficult not only to get access but also to measure its turgor pressure without damage. The first measurement of turgor pressure in the phloem was by Buttery and Boatman (Buttery and Boatman, 1964) using adapted manometers directly inserted in the laticiferous tissue of Pará rubber tree (Hevea brasiliensis). Although not belonging to the translocation pathway, the latex vessel system and sieve tubes are intimately associated elements so that the latex vessel system is considered part of the phloem tissue in Hevea and its turgor pressure was assumed to represent values in sieve tubes (Nicole et al., 1991). However, the first measurements of sieve tube turgor pressure (Hammel, 1968) used a manometric device, improved from that of Buttery and Boatman (1964), inserted in the secondary phloem of red Oak trunk (Quercus rubra). Hammel (1968) also observed gradients in osmotic pressure, pressure reducing at 0.03 MPa m −1 in the direction of flow from source to sink. However, those values were surprisingly low. According to the Münch hypothesis, higher values of turgor and osmotic pressure gradients were expected for the flows and sieve element anatomy of oak. Later, using pressure transducers and the phloem needle technique of Hammel (1968), other species were examined: white ash trunk (Fraxinus americana L.) (Sovonick-Dunford et al., 1981;Lee, 1981a,b); Manna ash (Fraxinus ornus) (Milburn, 1980); squirting cucumber stems (Ecballium elaterium) (Sheikholeslam and Currier, 1977a,b). By using a pressure bomb and Bourdon-type gauge Milburn and Zimmermann (1977) measured pressures in coconut palm inflorescences (Cocos nucifera L.). Milburn (1980), using a glass-spiral pressure gauge, measured pressure in castor bean (Ricinus communis) sieve tubes. Wright and Fisher (1980) developed another direct way of measuring sieve tubes turgor pressure using severed aphid stylets with attached capillary micromanometers to indicate pressure in sieve tubes of Babylon willow (Salix babylonica) stems and in bark strips of sandbar willow (Salix exigua) Wright and Fisher (1983). Following Fisher (1980), Fisher andCash-Clark (2000) measured sieve tube pressure in the peduncle of wheat. Following the same method, but with an on-line sensor, Gould et al. (2004) presented the second work on a smaller species sow thistle (Sonchus oleraceus). In an indirect method, phloem turgor pressure was determined by Kaufman and Kramer (1967) from measurements of tissue water potential along with phloem sap osmotic pressure in red maple (Acer rubrum L.), in reasonable agreement with Hammel (1968). Using the same method, Rogers and Peel (1975) observed turgor pressure gradients in stems of osier (Salix viminalis), and purple willow (Salix purpurea). Based on these findings, for our modeling purposes turgor pressure at entry to the system was set to p i = 1 MPa.

SIEVE TUBE MEMBRANE HYDRAULIC CONDUCTIVITY, L P
One of the main parameters describing sieve tube membrane transport is the membrane hydraulic conductivity, L p . There are not many data on L p . The first measurements, by Milburn (1974) on castor bean (Ricinus communis L.) bark segments, gave values between 5.7 and 8.8 × 10 −14 m s −1 Pa −1 . Sovonick-Dunford et al. (1982) obtained a value of 9.6 × 10 −15 m s −1 Pa −1 for the sieve tube elements of secondary phloem of red oak stem (Quercus rubra). Wright and Fisher (1983) measured 5 × 10 −15 m s −1 Pa −1 on sandbar willow (Salix exigua) bark strips. Before these investigations, most phloem transport models used the hydraulic conductivity published for membranes of other types of plant cells. However, most of that work was for algal cells. Tyree's (1970) review found the membrane hydraulic conductivity for plant cells to range from 5×10 −15 to 1 × 10 −12 m s −1 Pa −1 . Thus, when compared with other cells, it seems that sieve tubes are at the bottom of that range. In this study, sieve tube membrane hydraulic conductivity, L p , will be taken as 5 × 10 −15 m s −1 Pa −1 .

SOLUTE PERMEABILITY, P S , AND REFLECTION COEFFICIENT, σ
Of all the biological parameters used to describe phloem transport, the reflection coefficient of sieve tube membrane to solutes, σ, and the membrane permeability to those same solutes, P s , are the most difficult to determine. From an irreversiblewww.frontiersin.org December 2013 | Volume 4 | Article 531 | 13 thermodynamic analysis of solute and water transport across a membrane, Kedem and Katchalsky (1958) showed that the membrane parameters P s and σ should be related: V S is the specific volume of the solute. To our knowledge there are no measurements of these parameters for the sieve tube element/companion cell complexes. Hence, one can only speculate from studies of other types of plant cells. Diamond and Wright (1969) presented an empirical relation between the parameters for non-electrolytes from experimental data of Nitella mucromata cells. Although they do not present any theoretical justification, the data show reasonably well that σ decreases closely in parallel with increasing permeability P s . The permeability was smaller than 10 −8 m s −1 for σ greater than 0.8, and of the order of 10 −5 m s −1 for σ ∼ = 0. From a survey of experimental data on plant cells, Kargol et al. (1997) found a very good agreement with Equation (A.2) in measurements from maize roots, but not with measurements from isolated cells of Nitella translucens. Kargol et al. (1997) observed that the linear relationship of Kedem and Katchalsky agrees more with experimental data from artificial membranes. Kargol and Kargol (2000) showed a better agreement with experimental data when the right hand side of Equation (A.2) is multiplied by K = 0.021 ± 0.015. Subsequently, Kargol (2001) showed that K depends on the solute concentration within the membrane: K =cV S . However, as a first approach we will adopt the linear relationship (A.2) multiplying the right hand side by K = 0.031, which is closer to the average value obtained for the very few plant species presented by Kargol and Kargol (2000), however none of them referring to phloem tissue specifically. Inserting physiological conditions in sieve tubes (Table 1) into Equation (A.2), the solute permeability of the sieve tube membrane is given by: Hence P s varies between 10 −9 m s −1 , for totally permeable membrane (σ = 0), and 10 −11 m s −1 , for an almost impermeable membrane (σ = 0.99). Connor et al. (1977) and Legge (1985) measured vertical water potential gradients in mountain ash (Eucalytus regnans F. Muell.), between −0.007 and −0.034 MPa m −1 , with increasing height. This agrees with the estimated range of −0.01 to −0.03 MPa m −1 of the total water potential gradient necessary to drive transpiration stream in trees (Zimmermann and Brown, 1971). Regarding smaller plants, Begg and Turner (1970) measured apoplastic pressure gradients of 0.08 MPa m −1 in tobacco. There are not many studies on solute concentrations in the stem apoplast.    Connor et al. (1977) and Zimmermann and Brown (1971), we specified an apoplastic water potential gradient, d out dz , of 0.03 MPa m −1 and apoplastic osmotic pressure gradient, d out dz , 0.01 MPa m −1 , in our model. Therefore the apoplastic pressure gradient, dp out dz , surrounding the sieve tube will be 0.04 MPa m −1 , taking the direction of flow from apical to basal ends as is convenient for phloem flow.

APPENDIX B. MODEL BASED ON LOCAL APPLICATION OF HAGEN-POISSEUILLE EQUATION
Here we show that our zeroth order model Equations (27) and (31) can also be obtained using an approach similar to Thompson and Holbrook (2003a), based on conservation equations and a local application of the Hagen-Poiseuille equation together with appropriate boundary conditions. The same assumptions hold as described in Materials and Methods except for absence of a radial coordinate, i.e., the radial flow profile inside the tube is neglected because the solution inside the tube is regarded as well stirred in the radial direction. Also the diffusion of solutes within sieve tube sap is neglected.

VOLUME CONSERVATION
We start with considering an infinitesimal volume element [z − z, z] of the tube along the axial direction. The balance between incoming (axial) and outgoing (axial and lateral) fluxes in the infinitesimal volume element is given by J (z − z) = J (z) + J lat (z). Here J denotes axial volume flux [m 3 s −1 ] which is related to axial flux density j [m s −1 ] by J = A j, A = πR 2 being the cross section through which axial flux occurs. Lateral area of the volume element is A lat = 2πR z and lateral volume flux is J lat = A lat j lat . Thus In analogy to Equation (8) lateral flux density across the membrane is given by Starling's Equation (Nobel, 2009)  where the van't Hoff relation for dilute solutions has been inserted. Performing the limit z → 0, Equation (B.1) becomes

CONVECTIVE FLOW
Convective flow inside the tube is assumed to emanate from a pressure gradient according to a local applicability of the Hagen-Poiseuille equation such that axial flux density is given by with μ * = μ/β as before.

SOLUTE CONSERVATION
In steady state, incoming and outgoing solute fluxes for the small volume element of the tube have to be balanced: where the first two terms are axial solute flux into and out of the volume element, respectively. The other two terms describe lateral solute flux as sum of convective flux (solute dragged by solvent) and diffusive flux across the membrane. The concentration difference across the membrane has been abbreviated as c (z) = c (z) − c out (z) andc(z) denotes solute concentration in the membrane, a suitably defined function involving inside and outside concentration, e.g.

APPENDIX C. PREPARING MODEL EQUATIONS FOR NUMERICAL EVALUATION
In order to numerically evaluate the system of differential equations (27), (28) and (31)