1. Cell Biology
  2. Plant Biology
Download icon

From plasmodesma geometry to effective symplasmic permeability through biophysical modelling

  1. Eva E Deinum  Is a corresponding author
  2. Bela M Mulder
  3. Yoselin Benitez-Alfonso
  1. Wageningen University, Netherlands
  2. Institute AMOLF, Netherlands
  3. University of Leeds, United Kingdom
Research Article
  • Cited 2
  • Views 1,239
  • Annotations
Cite this article as: eLife 2019;8:e49000 doi: 10.7554/eLife.49000

Abstract

Regulation of molecular transport via intercellular channels called plasmodesmata (PDs) is important for both coordinating developmental and environmental responses among neighbouring cells, and isolating (groups of) cells to execute distinct programs. Cell-to-cell mobility of fluorescent molecules and PD dimensions (measured from electron micrographs) are both used as methods to predict PD transport capacity (i.e., effective symplasmic permeability), but often yield very different values. Here, we build a theoretical bridge between both experimental approaches by calculating the effective symplasmic permeability from a geometrical description of individual PDs and considering the flow towards them. We find that a dilated central region has the strongest impact in thick cell walls and that clustering of PDs into pit fields strongly reduces predicted permeabilities. Moreover, our open source multi-level model allows to predict PD dimensions matching measured permeabilities and add a functional interpretation to structural differences observed between PDs in different cell walls.

Introduction

The formation of spatial patterns in plants requires the transport and interaction of molecular signals. This sharing of information coordinates cell fate decisions over multiple cells and the isolation of cell fate determinants within a cell or group of cells on the same developmental path. Small molecules such as sugars, peptides, hormones and RNAs move long and short distances to coordinate cell/organ development (Otero et al., 2016). Cell-to-cell transport of proteins, such as transcription factors, is also important in the regulation and/or developmental reprogramming of local cellular domains (Gallagher et al., 2014). A well studied example is SHORT-ROOT (SHR), an Arabidopsis thaliana GRAS family transcription factor, that moves from the stele to cortical-endodermal tissue layers to specify cell fate and root patterning (Nakajima et al., 2001; Spiegelman et al., 2018; Wu and Gallagher, 2013; Wu and Gallagher, 2014). Other mobile factors with developmental importance include TARGET OF MONOPTEROS 7, PEAR transcription factors and miRNAs (Lu et al., 2018; Miyashima et al., 2019; Skopelitis et al., 2018).

Plant cells are connected by channels named plasmodesmata (PDs) that facilitate the transport of these molecules. PD are narrow membrane lined structures embedded in cell walls to allow for symplasmic (cytoplasm-to-cytoplasm) molecular flux (Figure 1). The ER forms a tubular structure called desmotubule (DT) that traverses the PD, leaving a discrete cytosolic sleeve (also called ‘cytoplasmic sleeve’ in the literature) where molecular transport occurs (Nicolas et al., 2017a; Sager and Lee, 2018). In the region closest to the PD entrances, the cytosolic sleeve appears constricted (neck) in most tissue types, although there are recent observations of ’straight’ PDs in meristematic root sections (Nicolas et al., 2017b). Cell walls at PD locations play a key role in regulating its dimensions. The accumulation of callose, a cell wall beta-1,3 glucan polysaccharide synthesized by callose synthases and degraded by β−1,3-glucanases (Zavaliev et al., 2011; Amsbury et al., 2017), is the best understood mechanism for the control of PD dimensions and symplasmic transport capacity (i.e. effective symplasmic permeability). Other factors such as membrane composition, shape and number of PDs change during development and between cell types adding extra dimensions to PD regulation (Nicolas et al., 2017a). Mutants blocked in PD form and function are embryo or seedling lethal, highlighting the importance of these structures for normal plant development (Kim et al., 2002; Benitez-Alfonso et al., 2009; Xu et al., 2012).

Modelling effective symplasmic permeability: concept overview.

(A) Electron microscopy image showing one PD, constricted at the neck regions (arrows), from Arabidopsis thaliana root tissue. The image was extracted from a reconstructed tomograph. Scale bar: 50 nm. The image was kindly provided by the Bayer lab. (B) Cartoon showing PD geometry and structural features. (C-F) The model to determine effective symplasmic permeability considers that connectivity within a cell file (C) is affected by the distribution of PDs in the cell wall (D) (modelled as a function of the cytoplasmic column belonging to a single PD (E)) as well as by the structural features of individual PDs (F).

Small molecules can move via PD by diffusion (non-targeted transport). This is considered to be predominantly symmetrical (Schönknecht et al., 2008; Maule, 2008), while in certain tissues, such as secreting trichomes (Waigmann and Zambryski, 1995; Gunning and Hughes, 1976) and the phloem (Ross-Elliott et al., 2017; Comtet et al., 2017), hydrodynamic flow may create directionality. The maximum size of molecules that can move by this generic ‘passive’ pathway is often referred to as the ‘size exclusion limit’ (SEL), which obviously depends on PD properties and structural features (Dashevskaya et al., 2008). Large molecules can move through PD via an 'active’ or ‘targeted’ pathway overriding the defined SEL. This may involve additional factors that temporarily modify these substrates, target them to the PDs, or induce transient modifications of the PDs to allow for the passage of larger molecules in a highly substrate dependent fashion (Zambryski and Crawford, 2000; Maule et al., 2011).

Computational modelling approaches have been applied to model PD transport but, so far, these have mainly focused on hydrodynamic flow and the specific tissues where that matters (Blake, 1978; Bret-Harte and Silk, 1994; Jensen et al., 2012; Ross-Elliott et al., 2017; Comtet et al., 2017; Foster and Miklavcic, 2017; Couvreur et al., 2018). The few existing studies on diffusive transport do not consider neck constrictions or the approach to PDs from the cytoplasmic bulk. Most models consider PDs as straight channels, with advective/diffusive transport through an unobstructed cytosolic sleeve and typically, but not always, account for reduced diffusivity inside these narrow channels compared to the cytosol (Bret-Harte and Silk, 1994; Liesche and Schulz, 2013; Dölger et al., 2014; Ross-Elliott et al., 2017; Couvreur et al., 2018). Only the oldest of this set, (Blake, 1978), uses a dilated central region in its calculations, but is entirely focused on hydrodynamics. In specific contexts, also a few other geometries are considered. (Ross-Elliott et al., 2017) also consider ‘funnel’ shaped PDs, which are observed in the phloem unloading zone, but ignore the DT in their diffusion model, as they only calculate a best case scenario for diffusive transport. In the context of size selectivity for small (sugar) molecules in phloem loading, also the so-called ‘sub-nano channel model’ of PD geometry has been considered (Liesche and Schulz, 2013; Comtet et al., 2017). In this model, symplasmic transport is modelled to be confined to nine cylindrical channels spanning the PD. This was based on a 9-fold rotational symmetry in enhanced 'top view’ electron micrographs but never validated experimentally in longitudinal sections. Instead, sparsely spaced axial spoke structures have been reported (Ding et al., 1992; Nicolas et al., 2017b).

Experimental measurement of the parameters that determine effective symplasmic permeability is difficult and many examples exist of misleading and/or conflicting results. Generally speaking two main approaches are used, providing results at different scales that are hard to reconcile. On the one hand, ultrastructural observations using transmission electron microscopy (EM) can provide useful data on PD dimensions and structural features but, despite recent advances, sample preparation affects the integrity and dimensions of PDs to an unknown extent potentially yielding an underestimation of relevant parameters (Nicolas et al., 2017b). On the other hand, tissue level measurement of symplasmic fluxes is achieved using symplasmic molecular reporters, but this is either invasive or limited to few molecular sizes, mostly fluorescein and its chemical relatives (hydrodynamic radii of about 0.4 to 0.6 nm) and GFP derived fluorescent proteins (such as DRONPA-s (28 kDa), Dendra2 (26 kDa), (photoactive and non-photoactive) single GFP (27 kDa, hydrodynamic radius 2.45–2.82 nm) and its multiples [Calvert et al., 2007; Terry et al., 1995; Chudakov et al., 2007; Gerlitz et al., 2018; Kim et al., 2005; Rutschow et al., 2011]). In all cases, the tissue geometry and varying degrees of vacuoloarization can severely complicate the interpretation of the measurements in terms of effective wall permeability for symplasmic transport. Old data on symplasmic permeability use either microinjection or particle bombardment, which allow for a much wider size range of dyes/molecular reporters, but these techniques can produce cellular stress, which affects PD function (Liesche and Schulz, 2012). Even when using the same dye/fluorescent molecule and the same tissue, these approaches deliver much lower permeabilities than less invasive techniques, demonstrating that they are unreliable for estimating permeabilities in unperturbed plants (e.g. see Haywood et al., 2002, or compare Rutschow et al., 2011 and Goodwin et al., 1990). Less invasive methods involve transgenic lines expressing fluorescent proteins under cell-specific promoters (Roberts et al., 2001; Stadler et al., 2005a), which are very time consuming to generate, or photoactivation and photobleaching techniques (Rutschow et al., 2011; Gerlitz et al., 2018). These approaches have yielded valuable insights, but again, both are limited to few proteins/molecular sizes.

In summary, despite recent advances in the development of probes and techniques, effective symplasmic permeability is difficult to assess directly. The fast response of plants to wounding and other stresses, may render part of the ultrastructurally derived parameters less reliable than others, explaining the frequent observation of apparently incompatible results when modelling diffusive symplasmic transport from ultrastructural measurements. In a multi-species analysis correlating photobleaching and electron microscopy results, (Liesche et al., 2019) were unable to find a universal model for matching measurements at the ultrastructural and tissue levels for different interfaces along the phloem loading pathway, illustrating the need for better models. Ideally, we would be able to integrate the results of the experimental approaches at both levels in a model that considers their limitations in order to get more accurate estimates of effective symplasmic permeability and the underlying structural parameters. This brings us to our central question: what do we need to assume about PD size, number, structure, etc. to be able to reproduce tissue level measurements? Moreover, PD geometry changes during development (Roberts et al., 2001; Fitzgibbon et al., 2013), inspiring our second main question: how do distinct features of PD geometry influence transport properties?

Here, we describe the biophysical properties of diffusive symplasmic transport considering detailed PD structural features (such as the DT and the neck region) and the approach from the cytoplasmic molecular bulk towards PDs that are either evenly distributed or clustered into pit fields (Faulkner et al., 2008) (Figure 1). Inside our model PDs, the entire cytosolic sleeve is available for particle diffusion (‘unobstructed cytosolic sleeve model’). We investigate how neck/central region, wall thickness, the presence of a DT and PD clustering into pit fields affect transport characteristics for different particle sizes, adding a functional context to some puzzling recent experimental observations. We also apply our framework to compute effective permeabilities for carboxyfluorescein (CF), a fluorescent dye used routinely to measure changes in symplasmic permeability. Comparing calculated and experimentally measured values, we demonstrate that the relatively high effective CF permeabilities observed by Rutschow et al. (2011) can be explained by our model of diffusive non-targeted symplasmic transport and reveal the potential source of conflicts with ultrastructural measurements. We found that, in this context, our model performed better than the ‘sub-nano channel model’ (Liesche and Schulz, 2013) referred to above. Our calculations demonstrate that multi-scale modelling approaches can integrate results from PD structural dimensions and molecular fluxes and reveal conflicts on these determinations. We, therefore, recommend these should be applied systematically when defining effective symplasmic permeability for a particular tissue/molecule and/or biological context. To facilitate this, we share a python program for computing effective permeabilities from PD geometries as a community resource.

Results

Outline of the model

Our aim is to describe the symplasmic transport properties of a cell wall as an effective wall permeability, that is a single number that could be plugged into tissue/organ level models. For this, we split the transport into two parts: the movement through an individual channel representing a PD and the approach to this channel from the cytoplasmic bulk (Figure 1). This implicitly assumes a homogeneous cytosol. The basic geometrical terminology that we considered in our calculations is introduced in the cartoon PD shown in Figure 1B. An overview of all mathematical symbols is given in Appendix 1.

Obtaining good EM data of PD dimensions is notoriously hard. We therefore opted for a simple geometrical description that allows us to study the effects of PD neck, central region and desmotubule dimensions with as few parameters as possible (see Materials and methods). We modelled a single PD as a 3-part cylindrical channel (Figure 2A), with total length l, which would typically equal the local wall thickness. The ends of the channel were modelled by narrow cylinders representing the plasmodesmal ‘neck’ constriction. These have length ln and radius Rn. The central region has radius Rc. Over the whole length, the center of the channel is occupied by a ‘desmotubule’ (DT) modelled as a cylinder of radius Rdt. The part available for diffusive transport, the cytosolic sleeve, is the space between the outer cylinder wall and the DT.

Model PD geometry and hindrance effects.

(A) Individual PDs are modelled using multiple cylinders with a total length l, neck (inner) radius Rn and neck length ln, central region (inner) radius Rc and desmotubule (outer) radius Rdt. B,C: Illustration of the impact of steric hindrance and rescaled parameters. The gray areas of the longitudinal (B) and transverse (C) sections cannot be reached by the center of the particle with radius α (steric hindrance). For a concise description of the available volume and cross section area, we use the rescaled lengths l~n=ln+α, R~c=Rc-α, R~dt=Rdt+α and R~n=Rn-α. (C) The cross section area available for diffusion on a transverse section was named A~, which depends on the particle radius (α). A~ is the area of the white ring in each cross section. The maximum particle size α¯ is illustrated with a dashed circle. For a particle of size α=α¯, A~=0. (D) In practice, particles spend less time diffusing close to the wall than farther away from it (hydrodynamic hindrance). Consequently, the area close to wall contributes less to diffusive transport, as illustrated with purple gradients. These additional hindrance effects are accounted for in A~~.

We made the arguably simplest choice of modelling particles as (non-additive, i.e. not interacting among themselves) hard spheres with radius α. This is partially supported by previous research showing that the hydrodynamic radius is the main determinant of PD transport characteristics, leaving behind, among others, particle charge (Dashevskaya et al., 2008; Terry and Robards, 1987). We also assumed that PD walls are rigid, and hence are unable to deform to accommodate larger particles. These assumptions imply a boundary condition: the center of a particle cannot come closer to the wall than the particle’s radius α (Figure 2B,C). This so-called steric hindrance reduces the volume that is available for diffusion of the particle’s center in a size dependent manner. Moreover, the maximum particle radius that can pass the PD, α¯, is always well defined. In practice, a precise definition of the SEL in terms of molecule size/shape is hard to give, however, we can use α¯ to operationalize the SEL concept in a straightforward manner. To avoid confusion, however, we will consistently write α¯ when referring to our model.

We introduced rescaled geometrical parameters to account for the reduced available volume in a compact way: l~n=ln+α, R~c=Rc-α, R~dt=Rdt+α and R~n=Rn-α. With these, the available surface area (Figure 2C) is

(1) A~x(α)=π(R~x2-R~dt2),(2α<Rx-Rdt),

with x=n for the neck and x=c for the central region. In the typical situation that the neck is the narrowest part of the channel, the maximum particle radius that can pass is: α¯=(Rn-Rdt)/2.

Considering pure diffusion without particle turnover inside the PD, particle flux through the channel is described by Cxyzt=D2Cxyz, or in steady state: D2Cxyz=0, with Cxyz the position dependent particle concentration and D the particle’s diffusion constant inside the PD. Note that D strongly depends on particle size. Assuming a homogeneous distribution of particle flux over (the available part of) each channel cross section, we can treat diffusion through the channel as a simple 1D problem along the channel axis (for the impact of this assumption, see Appendix 2). Particle mass conservation, as dictated by the steady state diffusion equation, then gives that the local concentration gradient at position x, Cx, is inversely proportional to the available surface area Ax, so Cc=A~n/A~cCn. The total concentration difference over the PD, ΔC=Cl-C0 is accordingly distributed over the channel: ΔC=2ln~Cn+(l-2ln~)Cc. The steady state molar flow rate Q(α) through each channel is proportional to the entrance cross section: Q(α)=-DA~nCn. Solving these equations for Cn leads to:

(2) Q(α)=-DA~nA~c2l~nA~c+(l-2l~n)A~nΔC.

This result can be improved further by incorporating hydrodynamic interactions between particles and walls (Figure 2D). To that end we followed (Liesche and Schulz, 2013) in employing the so-called hindrance factors 0H(λ)<1, which are based on proper cross sectional averaging of particle positions over time, as described by Dechadilok and Deen (2006). Based on geometrical considerations, we used the factors for a slit-pore geometry (see Materials and methods). These factors depend on the relative particle size λ. In our case, λ=2α/(Rx-Rdt). In the neck region, λ=α/α¯. For the full expression and behaviour of H(λ), see Materials and methods.

As H(λ) already includes the effect of steric hindrance between wall and particle, we can adjust Equation 2 by replacing every instance of Ax~ with

(3) A~~x=H(2αRx-Rdt)Ax.

For completeness, we note that the simplification of a uniform particle flux along the channel axis is violated near the neck-central region transitions, resulting in an error of a few percent (see Materials and methods for further discussion). We now define the permeation constant of a single PD, Π(α), through the rule rule steady-state flow rate = permeation constant × concentration difference, yielding

(4) Π(α)Q(α)ΔC=DA~~nA~~c2l~nA~~c+(l-2l~n)A~~n.

We also defined τ as the corresponding estimate for the mean residence time (MRT) in the channel. Using a steady state mass balance argument this can be calculated as the number of particles in the channel divided by the number leaving (or entering) per unit of time (see Materials and methods for further description).

(5) τ(α)=0lCxA~~xdx/Q(α)

Having defined the permeation constant of a single channel, the effective symplasmic permeability of the wall as a whole (P(α), the quantity that can be estimated using tissue level measurements) follows from the definition J=PΔC (steady state flux=permeability×density jump):

(6) P(α)=fihρΠ(α),

with ρ, the density of PDs per unit wall area (number/ μm2) and fih, a (density dependent) correction factor for the inhomogeneity of the wall (0<fih<1). The latter takes into account that the wall is, in fact, only permeable at discrete spots. To calculate fih, we considered a linear chain of cells of length L that are symplasmically connected over their transverse walls (Figure 1C) and computed mean first passage times (MFPT) through a straight PD and a column of cytoplasm surrounding the PD. The column was determined by assigning every bit of cytoplasm to the PD closest to it. For a regular triangular PD distribution, this results in a hexagonal column from the middle of one cell to the middle of the next, with a PD in its centre (Figure 1D). We then converted the MFPT to an effective wall permeability and compared the result with the uncorrected effective permeability computed as ρΠ(α) (as described in the Materials and methods).

As expected, P(α) depends on particle size. Two factors underlie this size dependence, which both affect Π(α): hindrance effects, which reduce the space available for particle diffusion, and the fact that the diffusion constant is inversely proportional to particle size: D=d1/α. Figure 3A and (Figure 3—figure supplement 1) show that hindrance effects have the strongest impact for particle sizes close to the maximum α¯, whereas the particle diffusion constant always has a large impact Figure 3B. For example, at Rn=Rc, the 50+ fold difference between α = 0.1 nm and α = 2 nm is reduced to a 3-fold difference when ignoring the particle size dependence of the diffusion constant.

Figure 3 with 1 supplement see all
Impact of particle size (radius = α) on single pore effective permeability Π(α).

(A) Dependence of Π(α) on neck radius (Rn) and α (different line colours, see legend). The diffusion constant D is inversely proportional to particle size (D=d1/α). Dashed lines show Π(α) considering only steric hindrance, solid lines include all hindrance effects. B: Using the same diffusion constant for all particle sizes instead shows that, once particles can pass easily, the particle size dependence of Π(α) is largely due to the relation between particle size and diffusion constant. Parameters for calculations: l = 200 nm, ln = 25nm, Rdt = 8 nm, Rc = 17.5 nm. For simplicity we use d1= 1 nm3/s in this figure. Therefore, only the relative values of the unit permeabilities have meaning (consequently expressed in arbitrary units [a.u.]).

Using the model presented here, we computed the effects of different PD structural features and changes in PD density and distribution on effective symplasmic permeability and its dependence on particle size as described below.

A dilated central region increases molecular flux in thicker cell walls

Electron microscopy suggests that PDs often have a neck region of reduced radius in comparison to the central region. We investigated how a constricted neck region, or, similarly, a dilated central region, affects PD transport. For this, we compared transport properties while conserving the size selectivity (constant α¯). We investigated how both the transport volume (using Equation 2) and transport time (τ as above) change when the central region is dilated. To compare channels with neck and dilated central region (12 nm =RnRc) with narrow straight channels (Rn=Rc= 12 nm), we define a relative molar flow rate as Qrel=Qdilated/Qnarrow and similarly relative τrel (Figure 4). For a more detailed discussion of τrel and its computation, see Materials and methods and Appendix 2.

Figure 4 with 2 supplements see all
Impact of central region dilation on molar flow rate (Q) and mean residence time (τ).

The same legend shown in C applies to all panels. Narrow channels have Rn=Rc= 12 nm, whereas for necked/dilated channels, Rn = 12 nm but Rc varies. (A-C) Red curves show the relation between molar flow rate in dilated PD vs narrow PD Qrel=Qdilated(Rn,Rc)/Qnarrow(Rn) whereas cyan curves show the relation between mean residence time in dilated PD vs narrow PD: τrel=τdilated(Rn,Rc)/τnarrow(Rn). Both quantities are computed for different particle sizes (solid: α0, dashed: α = 0.5 nm, sparse dashed: α = 1 nm, dash-dotted: α = 1.5 nm). (A, B) Qrel and τrel are shown as a function of the radius in the central region Rc for different PD lengths (cell wall thickness) (A) l = 100 nm, (B) l = 200 nm. (C) Values calculated for Rc = 17.5 nm (Rc* in A,B) as a function of PD length. (D) Ratios of curves calculated for Rc = 17.5 nm (C) and Rc = 26.4 nm (Figure 4—figure supplement 1B) represented for varying PD lengths. Other parameters used for modelling are: ln = 25 nm, Rn = 12 nm, Rdt = 8 nm.

We then investigated how both Qrel and τrel change with increasing central region radius Rc and how this depends on particle radius α and PD length l (Figure 4). The panels A and B show that molar flow rate increases with the central radius but quickly saturates, whereas mean resident time increases without upper bound. Moreover, both quantities increase faster for larger particle sizes (α, dashed lines). In fact, from studying the limiting behaviour of the underlying formulas, we found that Qrel is always less than its theoretical maximum l2l~n, whereas τrel ultimately scales quadratically with Rc, and, equivalently, linearly with the surface ratio A~~c/A~~n (see Appendix 3 and Appendix 3—figure 1). In simpler terms: the benefits of increased transport volume with increasing Rc saturate, and instead the costs in transport time increases ever faster with further dilation of the central region. This defines a trade-off between transport volume and transport time with increasing Rc when we analyze a single PD with a given entrance area.

Our computations also show that with increasing PD length l, the balance between both factors shift, because a much larger increase of Qrel is possible (Figure 4A–C). Similarly, for any given combination of Rn and Rc, Qrel decreases with increasing ln and decreases faster for shorter l, whereas τrel has its maximum at l~n=l/2 (Figure 4—figure supplement 2). Together, these computations suggests that dilation of the central region is more favourable in thicker cell walls. Interestingly, this theoretical observation correlates well with a recent EM study in Arabidopsis root tips (Nicolas et al., 2017b). The authors observed that PDs with a distinct dilated central region and neck region occurred mostly in thicker cell walls (average 200 nm), whereas in thin cell walls (average 100 nm), they found mostly straight PDs.

Additionally, (Nicolas et al., 2017b) observed a smaller and less variable radius in channels where the central region was occupied by spokes compared to channels without them (Rc = 17.6 nm vs. 26.4 nm on average). To analyze the effects of these changes on molar flow rate and MRT, we redrew the curves to compute relative values for Rc = 26.4 nm and Rc = 17.5 nm as a function of PD length (cell wall thickness) and for various particle sizes. As an example, panel C shows the variations observed when considering Rc = 17.5 nm (Rc* in A,B). We found that the molar flow rate Qrel increases less than the MRT τrel when increasing Rc from 17.5 nm to 26.4 nm, except for the smallest particle sizes in combination with large l (Figure 4D). These data suggest that in cell walls of moderate thickness, restricting the radius of the central region (which can be achieved by adding spokes) improves overall performance.

In summary, transport time and transport volume scale differently with the radius of the central region thus producing PDs with a dilated central region becomes more favourable when cell wall thickness increases. However, if the radius of the central region becomes too wide (as exemplified here for Rc = 26.4 nm) the increase in transport volume does not compensate for the delay in transport time. Interpretation of this result might explain why mostly straight PDs are found in recently divided cells (with thin cell walls) and why spokes (potentially limiting the radius of the central region) are often observed in mature PDs.

For the same given maximum particle size a PD with desmotubule can transport more than a PD without

A conserved feature of PDs –at least in embryophytes– is the presence of the DT, so we asked how this structure affects the transport capacity for particles of various sizes. In our model, the DT and the neck radius jointly define the maximum particle radius α¯. Assuming that control over maximum particle size α¯ is important and a high net flux often is desirable, we estimated the number of cylindrical channels required to match a single PD with DT. Using that P(α) is proportional to orifice area (An), we first computed nc(α¯), the number of circular channels that would offer the same An as a single channel with a DT of radius Rdt = 8 nm and the same α¯:

(7) nc(α¯)=(Rdt+2α¯)2-Rdt2α¯2=4Rdt+α¯α¯.

Figure 5A displays the nc(α¯) as a function of the maximum particle size. As an example, when α¯ = 2 nm, 20 cylindrical channels without DT would be needed to match the orifice surface area of a single channel with DT (with Rdt = 8 nm). This number decreases for larger α¯. We then considered that not all of this surface area is available for transport because of hindrance effects (Figure 2B–D). We found that even if the total surface area is the same, the channel with DT has a larger available surface area than the equivalent number of cylindrical channels. This is because in cylinders a larger fraction of the surface is close to the wall and, hence, hindrance effects are much more severe (Figure 5B, Figure 5—figure supplement 1). The difference increases with increasing relative particle size (α/α¯). Steric hindrance, that is the center of a hard particle cannot come closer to the wall than its own radius, plays only a minor part in this effect (Figure 5B).

Figure 5 with 1 supplement see all
DT increases the cross section surface area available for transport per channel given a maximum particle radius α¯.

(A) The number of cylindrical channels (nc) that is required to match the total entrance surface of a single channel with Rdt = 8 nm and the same maximum particle radius α¯. (B) Shows the relative area available for transport (An) in relation to relative particle size (α/α¯) when comparing channels with DT and the equivalent number of cylindrical channels. Total surface area is the same. Solid lines include all hindrance effects (A~~n/,dtA~~n,circle; cf. Figure 2D). Dashed lines includes steric effects only (A~n,dt/A~n,circle; cf. Figure 2C).

Clustering of PDs in pit fields reduces effective symplasmic permeability

The cell wall is only permeable for symplasmic transport where the PDs are. In this scenario, particles have to diffuse longer distances (on average) to reach a spot to cross the wall compared to a wall that is permeable everywhere. To account for this, we have introduced a correction factor, or ‘inhomogeneity factor’, fih in Equation 6 for the effective symplasmic permeability. Here, we explore how fih depends on all model parameters. To calculate fih, we treated the cytoplasm as a homogeneous medium. This simplifying assumption is necessitated by the lack of detailed information on the cytoplasm structure and how it differs among cells. Effectively, we assumed that the obstructing effects of ER, vacuoles, etc. are similar throughout the whole cell volume and thus can be captured in a single reduced cytoplasmic diffusion constant.

First, we calculated fih for isolated PDs positioned on a triangular grid in the cell wall (Figure 6A), as described in the Materials and methods. In Figure 6 we presented fih as a function of Rn and explored its dependence on particle size α (Figure 6—figure supplement 1A), presence/absence of DT (Figure 6—figure supplement 1A), cell length L (Figure 6—figure supplement 1B), density of PD ρ (B), wall thickness l (C) and PD distribution in the wall (D). We found that, provided that Rn is large enough for particles to enter (as indicated by vertical cyan lines in Figure 6—figure supplement 1A), fih is independent of cell length L and particle size α (Figure 6—figure supplement 1A,B) and is not affected by the DT. We also adjusted the computation for different regular trap distributions (Berezhkovskii et al., 2006) to find that fih also hardly depends on the precise layout of PDs (Figure 6D). Although variations in fih appear larger at low PD densities, for typical Rn values (for example, 12 nm as in Figure 4) density only has a minor impact (Figure 6B). Finally, we found an increase of fih with increasing PD length l, saturating to its theoretical maximum of fih=1 in thick cell walls (l > 500 nm) (Figure 6C). This result reflects the increasing time required for passing the PD itself with increasing PD length and, hence, a decreasing relative importance of the cytoplasmic diffusion.

Figure 6 with 1 supplement see all
Correction factor fih for inhomogeneous wall permeability depends on PD distribution, cell wall thickness and neck radius.

(A) The cartoon shows the geometrical considerations and parameters used to model the diffusion towards PDs. Cell wall inhomogeneity is incorporated as a correction factor fih, 0<fih1, which measures the relative impact of cytoplasmic diffusion towards the locations of the PDs in the cell wall compared to reaching a wall that is weakly but homogeneously permeable (i.e., with fih=1). The cytoplasm is considered homogeneous. Each bit of cytoplasm can be assigned to the PD closest to it. With PDs on a regular triangular grid, the cytoplasm belonging to a single PD, with an outer (neck) radius Rn, is a hexagonal column with cross section area Aw and 1/2 of the cell length L on either side of the wall. (B-D) fih is represented as a function of Rn. The presence/absence of DT does not affect the values of fih (Figure 6—figure supplement 1A). In all cases, solid lines correspond to: l = 100 nm, L = 10 μm, α = 0.5 nm, a PD density of ρ= 10 PD/μm2, and PDs distributed on a triangular grid. Broken lines show the effects of changes in PD density ρ (B), PD length l (C) and PD distribution (D).

Second, we investigated the effect of PDs grouped in small clusters resembling pit fields (see Materials and methods). The average centre-to-centre distance between PDs in pit fields considerably varies across species, with reported/calculated distances between 60 and 180 nm (Terauchi et al., 2015; Schmitz and Kühn, 1982; Danila et al., 2016; Faulkner et al., 2008). The lowest values, however, are from brown algae, which have a different PD structure from higher plants (Terauchi et al., 2012). As a default, we used d = 120 nm, which also coincides with measurements on electron micrographs of tobacco trichomes presented in Faulkner et al. (2008). In Figure 7A we calculated fih as a function of total PDs (‘entrances’) per area of cell wall for different numbers of PDs p clustered in a single pit field. We found that fih decreases with increasing number of PDs in a pit (and constant total PD density ρ). Different from isolated PDs, Figure 7A also reveals that, when grouped in pit fields, there is a strong dependence of fih on total PD density (number of PD entrances per area of cell wall). This could be predicted from extrapolating Figure 6B for isolated PDs, where density dependence also increases with increasing PD radius, because cluster radii Rpit are much larger than the largest Rn used in Figure 6B. Figure 7B shows that clustering (in this case 7 PDs) increases the dependence of fih on PD length (compare solid and dashed lines of the same colour). Increasing the distance between PDs within the cluster (Figure 7C), also increases the dependence of fih on PD density. Also the arrangement of PDs in small model clusters affects the degree of dependence fih on ρ. In both cases, we observe the steepest dependency of fih on ρ for the clusters with the lowest within cluster PD density (pit fields with p = 5, 6 and 19: indicated with blue lines in Figure 7A; see also Table 1).

Impact of PD clustering into pit fields.

PD organization within pits is indicated with small cartoons in each graph. Pits themselves are distributed on a regular triangular grid. Within pit fields, the nearest neighbour distance between PDs d (120 nm by default) is independent of the number of PDs per pit field. (A-C) fih is represented as a function of total PD density ρ (the total number of PD entrances per unit of cell wall area) for: a varying number of PDs per cluster p (as indicated by line type, (A), for different PD length l (B, solid lines: isolated PDs, dash-dotted lines: 7 PDs per cluster, red colour indicates l: 100 nm, cyan for 200 nm, blue for 500 nm) and for different PD spacing within clusters (C, shown for clusters of 7 PDs with centre-to-centre distance d as indicated by line type and colour). Cluster sizes 5, 6, and 19 are indicated with blue lines for readability (A,D). For comparison, fih for non-clustered but randomly distributed PDs is also indicated. (D) The impact of increasing the number of PDs per cluster p on P(α) as a function of cluster density ρpits (the number of pit fields per unit of cell wall area). Lines show the fold increase of P(α) when increasing the number of PDs per cluster from one to the number indicated by the line type (same as in A). Lines are terminated where fih of clusters meets fih of isolated PDs at the same total PD density. Beyond that, calculation results are no longer reliable because clusters get too close and the impact of clustering on fih could be considered negligible. (A-D) Default parameters: l = 100 nm, d = 120 nm, Rn = 12 nm.

Table 1
Pit radius (Rpit) as a function of number of PDs per pit.

The third and fourth column show numerical values for d = 120 nm and Rn = 12.

PDs/pitRpitAPD/ApitRpit
1Rn112
2Rn+12d0.05672
3Rn+133d0.06581.3
4*Rn+122d0.06196.9
5*Rn+d0.041132
6Rn+233d0.038150.6
7Rn+d0.058132
12Rn+1313d0.071156.2
19Rn+2d0.043252
  1. *: All entries are based on PDs on a triangular grid within each pit, except for 4 and 5, where the PDs inside a pit are arranged on a square grid. Clusters (pitfields) are always arranged on a triangular grid.

It is hypothesized that PD clustering arises or increases in the process of increasing PD number post cytokinesis, possibly through (repeated) ‘twinning’ of existing PDs (Faulkner et al., 2008). We, therefore, also investigated the effect of increasing the number of PDs per cluster (p), starting from 1 PD per cluster (Figure 7D). As expected, P(α) always increased with the increase in cluster size/PD number (Figure 7D), despite the decrease in fih compared to homogeneously distributed PDs. This increase was larger for larger pit densities (number of pit fields per cell wall area).

In summary, for isolated and roughly evenly distributed PDs, the correction factor fih for inhomogeneous wall permeability has only a minor role on P(α). For realistic PD dimensions (Rn < 20–25 nm), the additional effect of fih with parameter changes would be too small to be observed experimentally, with the possible exception of PD length l. However, when considering clusters of PDs, as is common in pit fields, fih is markedly reduced, and PD length and density have a much larger impact on fih. We observed the biggest difference between isolated PDs and pairs, that is when going from single to twinned PDs (Figure 7A).

Application of the model to compute effective permeability for fluorescein derivatives

In a system where non-targeted symplasmic transport is fully driven by diffusion (so no (significant) active transport or hydrodynamic flow), our calculations using reasonable PD dimensions and densities should yield values close to the ones measured experimentally. As a resource to test this hypothesis, we have build a Python program, PDinsight, that computes effective permeabilities from parameter measurements extracted from EM. As some of these parameters might be more reliable than others, we also created a mode in the program to predict what are the minimum requirements in terms of parameter (combination of parameters) values to obtain experimentally measured symplasmic permeability. Exploring these requirements is equivalent to testing hypotheses like: ‘What if PD aperture is larger than observed with EM? or if the molecular radius is smaller than predicted?”. Predictions made with the program can be used to explain experimental results, highlight areas/parameters that need more investigation and can help with the design of new strategies to change effective symplasmic permeability in vivo. For a full description of the program and its possibilities, see Appendix 6.

As a test case, we used the model to explain the permeability measurements in Arabidopsis thaliana roots reported for carboxyfluorescein (CF) diacetate: a membrane permeable non-fluorescent dye that once converted inside cells into a fluorescent version of fluorescein can only move from cell to cell via the PDs (Rutschow et al., 2011). Using a technique named fluorescence recovery after photobleaching (FRAP), CF effective permeability was estimated for transverse walls in the root meristem zone (measured ≈ 200 μm from the quiescent centre). The authors present two experimental setups: a ‘tissue level’ experiment in which a whole ≈ 50 μm longitudinal section of the root was bleached (estimated effective permeability 6–8.5 μm/s) and a single cell experiment in which a single epidermal cell was bleached (estimated effective permeability 3.3 ± 0.8 μm/s).

PD densities in transverse walls of Arabidopsis thaliana roots were reported by Zhu et al. (1998): vascular: 9.92 ± 0.58, inner cortex: 12.28 ± 0.67, outer cortex: 9.08 ± 0.50 and epidermis 5.42 ± 0.42 PDs/μm2. Based on these numbers we assume a PD density of 10–13 PDs/μm2 for the tissue level experiment and 5 PDs/μmfor the single cell experiment. Fluorescein has a Stokes radius of approximately 0.5 nm (Champion et al., 1995; Corti et al., 2008) and a cytoplasmic diffusion constant of D = 162 μm2/s (one third of its water value) (Rutschow et al., 2011). Feeding these numbers to the model, and considering that PDs appear as straight channels in these walls (Nicolas et al., 2017b), we are able to reproduce the measured permeability values for observed PD densities (Zhu et al., 1998) only if we assume a relatively wide open neck (Rn > 15 nm) (Figure 8A,B, Table 2). Notably, the required neck radius for the single cell experiment fits within the range of the tissue level experiment when considering the respectively measured densities. This prediction is plausible if we consider that, in the same tissues, GFP (a protein with a reported hydrodynamic radius of 2.45 nm [Calvert et al., 2007] to 2.82 nm [Terry et al., 1995]) moves intercellularly (Stadler et al., 2005b). Using our default Rdt, Rn should be distinctly wider than 13–14 nm for GFP to move. We also explored the possibility that PD densities are higher than determined by Zhu et al. (1998). We found that to obtain the required effective permeabilities for CF with our default Rn = 12 nm, we would need PD densities of 33–47 PDs μm-2 for the tissue level experiment and 19 (14 - 23) PDs μm-2 for the single cell experiment (Table 2). The ratio of these required densities is in line with the observed ratio of relevant densities (Zhu et al., 1998).

Calculated effective permeabilities for carboxyfluorescein (CF) as a function of PD aperture at the neck Rn.

(A, B) Shows the graphs for straight channels. (A) Effective permeabilities are calculated for different PD densities (different colour curves). The horizontal gray band in A and C indicates the cortical values observed by Rutschow et al. (2011). (B) Shows the PD density required to obtain measured values of P(CF) (different colour curves) as a function of Rn. Horizontal broken lines are introduced to aid readability. (C, D) Shows that effective permeability increases with dilation of the central region (Rc>Rn). As a reference, values for straight channels are indicated in black. Dashed curves show values calculated for channels without DT. (D) Shows the same calculations as C but for longer PDs l = 200 nm. Default parameters: α = 0.5 nm, D = 162 μm2/s, ln = 25 nm, l = 100 nm, Rdt = 8 nm, ρ = 10 PD/μm2, PDs are spaced on a triangular grid, without clustering.

Table 2
Parameter requirements for reproducing measured P(CF) values (Rutschow et al., 2011) with the default model.

This table was generated using PDinsight. A: Required density (ρ) for a given α¯ and neck radius (Rn). B: Required α¯ and corresponding Rn for a given ρ. C: values required to reproduce P(CF) = 25 μm/s. Values computed for a 2x, 3x and 4x increase of ρ are also shown. This is done both for a uniform increase of the density (p=1) and for (repeated) twinning (p>1) from a uniform starting density (indicated in bold). p is the number of PDs per pit.

P(CF)(μm/s)α¯ (nm)Rn (nm)ρ (PD/μm2)
3.3*2.01218.6
2.51312.6
3.0149.3
3.414.87.6
3.5157.3
4.0165.9
62.01233.5
2.51322.7
3.01416.8
3.414.813.8
3.51513.2
4.01610.7
8.52.01247.2
2.51332.0
3.01423.7
3.414.819.4
3.51518.5
4.01615.0
P(CF)(μm/s)ρα¯Rn (nm)
3.3*54.516.9
4.216.5
6104.216.3
133.515.1
8.5105.218.4
134.416.8
1101.511.0
131.310.6
P(CF)(μm/s)ρpα¯Rn (nm)
2510110.528.9
2016.621.2
27.222.5
3015.118.1
35.619.2
4014.216.4
44.617.2
1318.825.6
2615.619.1
26.020.0
3914.316.6
34.617.2
5213.615.1
43.815.5
  1. *: Single cell experiment. All other data relates to tissue level experiments.

Using the model, we also explored the effect of ‘necked’/‘dilated’ PDs by adding a wider central region to PDs. For a central radius Rc = 17.6 nm, the required Rn to reproduce the tissue level CF permeability values would decrease by perhaps 1 nm or at most 3 nm (for Rc = 26.4 nm) considering a PD density in the order of ρ = 10 μm–2 (Figure 8C, Rc values from Nicolas et al., 2017b). In thicker cell walls (Figure 8D), the calculated effective permeabilities increased relatively more, but remained too low, suggesting that increasing cavity radius is never sufficient for reproducing the Rutschow et al. (2011) values (see also Figure 4).

Using the tissue level setup, Rutschow et al. also reported drastic changes in effective permeability after H2O2 treatment. They found a strong decrease in symplasmic permeability to ≈ 1μm/s after treatment with a ‘high’ H2O2 concentration, which was explained by rapid PD closure through callose deposition. Using our program we found that, for this reduction of P(CF), callose must reduce Rn to 11 nm (ρ = 10 μm–2) or 10.6 nm (ρ = 13 μm–2), resulting in α¯ = 1.5 nm or 1.3 nm, respectively. The authors also found a strong increase in permeability to ≈ 25 μm/s after treatment with a ‘low’ H2O2 concentration. Reproducing this increase requires a large change at the PD level. At the extremes, an increase of Rn to approximately 29 nm for ρ = 10 μm–2 (Figure 8A,B, Table 2B), or a slightly more than four fold increase in PD density would be required to reproduce this high effective permeability (Table 2C). Alternatively, both Rn and ρ would have to increase substantially (Figure 8B). As an extreme hypothesis, we also calculated the effects of complete DT removal. The increases in P(α) that could be obtained this way were by far insufficient to explain the reported effect of mild H2O2 treatment (Figure 8C,D), making DT modification or removal a highly unlikely explanation for this change.

Taken together, these calculations indicate that our model for diffusive symplasmic transport can indeed explain experimentally observed measurements of effective symplasmic permeability, but only with somewhat wider PDs/neck regions than expected yet in line with the observed permeability for GFP and within the range of PD diameters measured in thick cell walls. Alternatively, similar changes in symplasmic permeability can be achieved with several fold higher densities than typically measured. These predictions provide a framework for experimental validation. We also compared the results obtained with our unobstructed sleeve model and the sub-nano channel architecture. Using the sub-nanochannel architecture, much larger PD densities would be required to achieve the same P(CF): roughly twice as large for α¯ = 3.5–4 nm and even larger for smaller α¯ (see Appendix 5 and Appendix 5—table 1). These results favour unobstructed sleeve models for offering more plausible hypotheses to explain the experimental results for CF and the impact of H2O2 treatment on effective permeability.

Discussion

In this manuscript, we presented a method for calculating effective wall permeabilities for non-targeted, diffusive symplasmic transport based on the dimensions and distribution of PDs and on the size of the mobile particles. For individual PDs, we used a minimal geometrical description that allowed us to extensively investigate the effects of dilation of the central PD region and the implications of a DT at the PD axis on transport properties. Because PDs are narrow, our calculated effective symplasmic permeabilities were heavily affected by molecular hindrance effects. For the effects of PD distribution, we introduced an ‘inhomogeneity factor’ fih between 0 and 1, which accounts for the reduction in overall permeability due to spatial arrangement of PDs. We found that the degree of PD clustering had a strong impact on this factor, whereas the exact spatial distribution of either isolated PDs or clusters had little impact.

Our model uses an unobstructed cytosolic sleeve for symplasmic transport. In such models, the DT gives the PD an annular cross section, which strongly increases transport capacity compared to cylindrical channels with the same α¯ and total cross section area at the entrance, particularly for relatively large molecules. Having a DT offers an additional flexibility in regulating size selectivity through the possibility of a dilated state of the PD by displacement or temporal removal of the DT (Zambryski and Crawford, 2000; Crawford and Zambryski, 2000). This feature, however, can be exploited for the spreading of viruses (Benitez-Alfonso et al., 2010) and other intracellular parasites such as the fungus Magnaporthe oryzae (Kankanala et al., 2007). Functional PDs without DT (and inner diameter of 10–20 nm) have been reported for the brown algae species Dictyota dichotoma (Terauchi et al., 2012). Due to their very high membrane curvature, DT formation requires curvature-inducing proteins (such as reticulons) and a special lipid composition (Tilsner et al., 2011; Grison et al., 2015; Knox et al., 2015). It is likely that performance benefits of the DT offset these costs and disadvantages and it is therefore under evolutionary selection. Additionally, the connection between DT and ER could result in variable degrees of PD occlusion and hence a potential control mechanism for PD accessibility. Park et al. (2019) have started to explore this concept in the context of pressure regulated PD occlusion.

We have also calculated the performance costs (transport rate) and benefits (transport volume per PD) of having distinct central and neck regions. Whereas the transport time scales quadratically with the radius of the central region (Rc), the relative transport volume has a strong upper bound that increases with channel length. These results suggest that straight PDs perform better in thin (average 100 nm) cell walls and necked/dilated PDs in thick (average 200 nm or more) cell walls, which correlates with recent observations (Nicolas et al., 2017b). This is not, however, the only way to explain these observations. Necked/dilated PDs might appear because (1) size selectivity is more efficiently controlled by restricting callose deposition to a 20–30 nm long neck region, (2) the formation of ‘spokes’ in the central region leads to this narrow-wide-narrow structure, and/or (3) the material properties of cellulosic cell walls and PD cell membranes only allow for a distinctly wider central region if the channel is long enough.

In our model, we naturally define the SEL as α¯, the maximum particle radius that could fit through the model PD, but experimental determination of this value is difficult and often relies on the transport of detectable, typically fluorescent, molecules such as CF. The limited set of suitable molecules, particularly for non-invasive techniques, introduces a large uncertainty in SEL measurements and hence α¯. Also other biological factors could lead to an underestimation as well as an overestimation of α¯. For example, in so-called active symplasmic phloem loaders, such as the cucurbits, sucrose moves symplasmically from bundle sheet cells (BSC) to intermediary cells (IC), where it is polymerized into the larger oligomers raffinose and stachyose, that do not diffuse back in detectable amounts (Haritatos et al., 1996; Liesche and Patrick, 2017). Two explanations have been suggested: (1) a discriminating PD SEL at this interface, which prevents the back transport of raffinose and stachyose (Liesche and Schulz, 2013), or (2) open PDs combined with a directional flow which could be sustained by the xylem flow (Comtet et al., 2017). Only the latter could explain the observed amount of sucrose transport (Liesche and Schulz, 2013; Comtet et al., 2017). This example illustrates that the consideration of a symplasmic flow could largely affect calculated permeabilities and fluxes.

An overestimation of α¯ could occur for non-spherical molecules or temporal variations in PD properties. Although a molecule’s hydrodynamic radius is a better predictor of its symplasmic transport efficiency than its molecular weight (Terry and Robards, 1987; Dashevskaya et al., 2008), it conceptually assumes a static replacement sphere. Molecules may be more flexible and/or have a shortest dimension than what is captured by its diffusive behaviour in bulk. PDs might also accommodate molecules that are larger than expected, either through interactions with specific PD proteins (Benitez-Alfonso et al., 2010) or because membranes and/or cell wall domains around PDs allow for reversible transient modifications in α¯ (Abou-Saleh et al., 2018). Additionally, molecules could pass in the wake of larger proteins/complexes/structures that modify PDs (e.g., tubule-forming viruses; Amari et al., 2010). Assessing the extent and time scales of temporal variations in PD boundaries and their implications remains an open topic for future investigation.

The framework we have developed for so-called ‘simple’ PDs also provides an intuition for the functional implications of complex geometries such as ‘twinned’, ‘branched’ or ‘funnel’ PDs (Ehlers and Kollmann, 2001; Ehlers and van Bel, 2010; Faulkner et al., 2008; Ross-Elliott et al., 2017). All else remaining equal, ‘twinned’ PDs have twice the entrance surface area, which would result in doubling the effective permeability P(α). This increase, however, will be reduced because of the less uniform PD spacing in a density dependent manner (Figure 7A). ‘Branched’ or ‘complex’ PDs contain multiple sub-channels (branches) on at least one side with typically a single shared central cavity connecting all branches (Oparka et al., 1999; Roberts et al., 2001; Fitzgibbon et al., 2013). In the leaf sink/source transition, massive branching is observed and, coincidentally, the number of PDs is reduced (Roberts et al., 2001). The formation of many channels per PD could help to maintain sufficient transport capacity for smaller molecules. If so, the increase in the number of typically narrower channels should be much larger than the decrease in total (simple or complex) PD number. Our computations of fih after twinning suggest that minimizing the distance between sub-channels could be favourable at low to moderate PD densities (Figure 7C). ‘Funnel’ PDs are reported in tissues surrounding the phloem at the root unloading zone (Ross-Elliott et al., 2017) and show a wide opening on the PSE (protophloem sieve element) side and a narrow opening on the PPP (phloem pole pericycle) side. (Ross-Elliott et al., 2017) model these as a triangular funnel that reaches its narrowest diameter only at the (PPP) bottom. There appears to be, however, a longer neck-like region at the narrow end of variable length. As hindrance is by far the highest in the narrowest section, the length of this narrow part would be a vital parameter in correctly estimating the transport permeabilities of these PDs.

We have applied our model to calculate the effective permeability for fluorescein in transverse walls of Arabidopsis root tip cells (Rutschow et al., 2011). Assuming purely diffusive transport and parameters based on various ultrastructural measurements, we were able to reproduce the observed effective permeabilities for CF and to assess the plausibility of different hypotheses aimed at resolving the conundrum of apparently incompatible measurements at different scales. For resolving this conundrum, we assumed that not all PD dimensions are reliably measured with EM. We could reproduce the measured values with somewhat wider PDs/neck regions or several fold higher PD densities than usually measured by EM. Of these, the increased radius seems the more plausible scenario, in line with the requirements for efficient GFP transport reported to occur among root meristem cells (Benitez-Alfonso et al., 2009; Benitez-Alfonso et al., 2013; Nicolas et al., 2017b), and similar to Rc values reported in thicker cell walls (Zhu et al., 1998; Grison et al., 2015; Nicolas et al., 2017b). Remarkably, our model predicts very similar PD aperture in the transverse walls of the epidermis and the more interior root layers when considering the ≈ 2-fold difference in PD density (Zhu et al., 1998). The obvious next step would be testing more data sets of different interfaces/plant species where purely diffusive symplasmic transport is expected. First of all, it would be ideal to test if a near or complete match between tissue level and ultrastructural measurements can be produced if all measurements are performed on the same system with the same growth/treatment conditions. Additionally, more testing could yield a better understanding of potential systematic side-effects of modern EM preparation techniques and/or uncertainties in the tissue level measurements, which would show as systematic vs random required adjustments of the model parameters. A very exciting outcome would be the discovery of distinct clusters in required parameter adjustments that could be related to cell wall properties, PD or interface type, etc. Additional model testing would become easier if the results of tissue level experiments are reported in the form of effective symplasmic wall permeabilities (in μm/s), or clearly provide all information required to transform into such units.

We also used our model to predict the PD changes after treatment with high and low concentrations of H2O2 in Rutschow et al. (2011). The reduced permeability after high H2O2 treatment could easily be explained by a redox induced stress response and corresponding reduction of PD aperture (e.g., at a density of 10 PD/μm2, a reduction from α¯ = 4.2–5.2 nm to α¯ = 1.5 nm would be required, see Table 2B). The strongly increased permeability after low H2O2 treatment, however, is harder to explain. With a single parameter change, the model predicts either a very wide PD aperture of α¯ = 8.8–10.5 nm, or a ±4-fold increase in PD density (possibly through 2 rounds of twinning/duplication), or less extreme changes if both parameters increase simultaneously (see Table 2C). The required increase in PD density should occur relatively fast, that is within the applied incubation period of 2 hr, and is so large that it should be readily detectable with EM.

The fact is that to reproduce experimentally measured CF effective permeabilities with our model, we had to deviate from ultrastructural based values for at least one parameter. Potential sources for these variations are: (1) ultrastructural studies might underestimate Rn because plants could respond to pre-EM manipulation by closing PDs, like they do in response to microinjection or particle bombardment (Haywood et al., 2002; Liesche and Schulz, 2012), (2) PD integrity could be affected during processing for TEM leading to an underestimation of PD densities, (3) the mechanical properties of cell walls and membranes provide a flexibility in the channel that could to some degree accommodate molecules larger than the apparent Rn (Abou-Saleh et al., 2018; Yan et al., 2019; Amsbury and Benitez-Alfonso, 2019). For a passive transport mechanism, the elastic energy required for these reversible deformations would have to be in the order of a few kBT or less. A model with flexible PD lining would be required to investigate the physical limits of this ‘flexibility hypothesis’, which is quite an increase in model complexity compared to the hard walls used in all current models, including ours. Finally, technical issues limit the accuracy of the CF effective permeability measurements themselves, for example, the speed of confocal microscopy bounds the spatial and temporal resolution at which CF concentrations can be monitored during and after bleaching/photoactivation (Rutschow et al., 2011; Liesche and Schulz, 2012).

To assess the impact on effective symplasmic permeability of various PD distributions, including clustering into pit fields, we introduced the inhomogeneity factor fih that accounts for the fact that the wall is only permeable at certain spots (i.e., where the PDs are located). Clustering into pit fields had by far the largest impact on this factor, particularly for lower PD densities. This means that not only total PD density, but also the degree of clustering is important information for calculating effective wall permeability from experimental data. The above inhomogeneity factor and the possibility of a dilated central region set our model apart from other models based on the unobstructed sleeve architecture (Bret-Harte and Silk, 1994; Liesche and Schulz, 2013; Dölger et al., 2014; Ross-Elliott et al., 2017). Using typical PD dimensions and no clustering, inhomogeneity factor fih would reduce the effective symplasmic permeability by about 15%, meaning that our model would require slightly wider or more PDs to explain the same tissue level experiments with straight channels compared to the above models.

A dilated central region is also considered in Blake (1978), who investigates hydrodynamic flow only. There is, however, an interesting similarity between both conditions: in both cases the driving gradient is steepest in the (narrowest part of) the neck region, be it the concentration gradient (Appendix 2—figure 2A) or the pressure gradient (Blake, 1978). When it comes to describing the PD geometry, (Blake, 1978), makes the opposite choice compared to us. He glues together sin2 functions with a straight middle part, resulting in a mathematically nice (i.e., continuous differentiable) function, but consequently, neck shape cannot be controlled, and neck length and the length of the widening region are linked. We, on the other hand, use an instantaneous increase in PD radius, which introduces a mild systematic error in our estimates of effective symplasmic permeability P(α) (Appendix 2), but results in parameters that are directly measurable on EM images.

Comparing the unobstructed sleeve architecture to the sub-nano channel architecture, we found that the latter requires roughly twice as high PD densities to produce the same permeability values P(CF) in the (Rutschow et al., 2011) experiments. This difference is due to the increased hindrance effects in cylindrical channels vs annular channels with the same cross sectional area. In the future, sleeve models could be refined with the consideration of central spokes (Ding et al., 1992; Nicolas et al., 2017b) and variability of PD dimensions within a single cell wall (Nicolas et al., 2017b; Yan et al., 2019). Simple considerations of the available volume suggest that the addition of spokes will increase hindrance effects, but most likely to a lesser extend than the sub-nano channel structure. Detailed molecular simulations could be a valuable tool to assess this effect.

Other future applications could be the coupling of our detailed PD level calculations of effective symplasmic permeability with tissue level models, which would allow for investigating the impact of microscopic changes on developmental and physiological processes (for example see Foster and Miklavcic, 2017; Couvreur et al., 2018). Depending on the context, it would then be useful or even required to also implement hydrodynamic flow through the PDs. Many ingredients are available for doing this while maintaining the distinguishing features of our mode, including hindrance factors (Dechadilok and Deen, 2006), but as far as we know, the theoretical and numerical results that we use for calculating fih are only available for diffusion processes, and not yet for advection. Additionally, one may need to replace the abrupt change in PD radius by a more gradual function. The importance of this final change could be estimated using numerical simulations.

Technological advances have started to be applied for more refined determinations on ultrastructural parameters. New fixation and sectioning techniques and new technologies such as electron-tomography (ET) and Correlative Fluorescence Electron Microscopy (CFEM) are now part of the systematic study of PD connections in different plant cells, tissues and organs. In parallel, new information on structural features characterizing PDs in different plant species/developmental stages as well as on the factors controlling PD structure and function (and thereby the effective permeability of specific molecules in different developmental or environmental conditions) are emerging. Combined with this significant experimental progress, our calculations provide a functional interpretation to characteristic PD morphological features and provide a framework to investigate how transport properties depend on these ultrastructural features and particle size in the context of simple and complex PD geometries. Another level of predictive power could be unlocked by integrating our framework into larger models at the tissue to whole organism level. This opens new avenues for exploring how developmental regulation of symplasmic transport interacts with various other pathways for long and short range intercellular communication.

Materials and methods

Diffusive flux through a single PD

Request a detailed protocol

Similar to Smith (1986), we assumed the flux is distributed homogeneously within each cross section along the axis of the channel. This results in a simple mapping to a 1D channel, that is that the average local flux (per unit area of cross section) ∼ 1/available cross section surface. This assumption does not hold close to the transition between neck and central region, that is a sharp transition between narrow and wide cylinders. Numerical simulations showed, however, that the error introduced by the assumption of homogeneous flux turned out to less than 4 percent for l = 200 nm, the shortest l with experimentally observed neck region in Nicolas et al. (2017b) (Appendix 2—figure 1) and will be less for longer channels. This error can be considered irrelevant given the quality of available data on PD dimensions and the many molecular aspects of PD functioning that are necessarily neglected in a simple model.

Hindrance factors

Request a detailed protocol

Hindrance factors H(λ) including both steric and hydrodynamic effects are modelled using the numerical approximations in Dechadilok and Deen (2006). They present functions for cylindrical and slit pores. For PDs with a desmotubule, we use the function calculated for straight slits.

(8) H(λ)=1+916λln(λ)1.19358λ+0.4285λ30.3192λ4+0.08428λ5.

This choice is supported by the steric hindrance prefactor that is included in H(λ) (Dechadilok and Deen, 2006). This Φ(λ)=1-λ is the same as the ratio of available to full surface area A~x(α)/Ax. For cylindrical channels, that is reference channels in Figure 5 and the regular PDs after DT removal, we use

(9) Hc(λ)=1+98λlog(λ)1.56034λ+0.528155λ2+1.91521λ32.81903λ4+0.270788λ5+1.10115λ60.435933λ7

for λ<0.95 and the asymptotic approximation by Mavrovouniotis and Brenner (1988),

(10) Hc(λ)=(1-λ)2(0.984(1-λλ)52

otherwise, as suggested by Dechadilok and Deen (2006).

Relative molar flow rate and MRT

Request a detailed protocol

For assessing the impact of the neck constriction on PD transport, we defined two relative quantities: Qrel=Qdilated/Qnarrow and τrel=τdilated/τnarrow (Figure 4, Appendix 3—figure 1). Using Equation 2 for Q(α), Qrel is well defined:

(11) Qrel(α,Rc)=lA~~c2(l~n)A~~c+(l2l~n)A~~n
(12) limRcQrel(α,Rc)=l2l~n

For τrel we first needed an expression for τ itself. Ideally, this would be a MFPT, which could calculated in a way similar to τ in the calculation of fih, using a narrow-wide-narrow setup. These calculations, however, critically depend on trapping rates at the narrow-wide transitions. We do not have an expression for these, because the DT takes up the central space of the channel, which, contrary to the case of fih, substantially alters the problem and the circular trap based calculations would result in an underestimation of the MFPT. Instead, we stuck to the homogeneous flux assumption also used for Q(α) and defined τ as the corresponding estimate for the mean residence time (MRT) in the channel (see Equation 5). Elaborating Equation 5:

(13)τ(α)=Cl+C02DΔC(2l~nA~~n+(l2l~n)A~~c)(2l~nA~~c+(l2l~n)A~~n)A~~nA~~c(14)=Cl+C02DΔC(4l~n2+(l2l~n)2+2l~n(l2l~n)(A~~cA~~n+A~~nA~~c)).

Unfortunately, this depends on the concentration difference over the channel. We are interested, however, in how the MRT changes with increasing Rc. In our definition of τrel, the concentration difference cancels from the equation, solving the problem:

(15) τrel(α,Rc)=1l2(4l~n2+(l2l~n)2+2l~n(l2l~n)(A~~cA~~n+A~~nA~~c)).

This method of computing τrel again depends on the homogeneous flux assumption. For an estimate of the error introduced by this approach, see Appendix 2.

Flow towards PDs: correction for inhomogeneity of the wall permeability

Request a detailed protocol

To compute fih, we consider a linear chain of cells that are symplasmically connected over their transverse walls (Figure 1). We first compute mean first passage time (MFPT) τ through a simplified PD and a column of cytoplasm surrounding it. We then convert τ to an effective wall permeability and compare the result with the uncorrected effected permeability computed using Equation 6 for the simplified PD geometry and fih=1.

As a simplified PD, we use a narrow cylindrical channel of length l and radius Rn, that is initially without DT. We assume that PDs are regularly spaced on a triangular grid. Consequently, the domain of cytoplasm belonging to each PD is a hexagonal column of length L, the length of the cell (Figure 6). We adjust the results reported by Makhnovskii et al. (2010) for cylindrical tubes with alternating diameter by changing the wide cylinder of radius Rw with a hexagonal column with cross section area Aw=1/ρ and considering hindrance effects. Makhnovskii et al. use a setup with an absorbing plane in the middle of a wide section and a reflecting plane, where also the initial source is located, in the middle of the next wide section. Assuming equal diffusion constants in both sections, they report the following MFPT from plane to plane:

(16) τ=12D[L2+l2+2D(lκn+Lκw)+lL(κwκn+κnκw)],

where

(17) κw=4DRnf(Rn2Rw2)πRw2

is a trapping rate to map the 3D setup onto a 1D diffusion problem. In this,

(18) f(σ)=1+Aσ-Bσ2(1-σ)2

is a function that monotonically increases from 0 to infinity as σ, the fraction of the wall occupied by the circular PDs, increases from 0 to 1. f(σ) is the result of a computer assisted boundary homogenization procedure with the values of A and B depending on the arrangement of trapping patches (Berezhkovskii et al., 2006). To maintain detailed balance, the corresponding trapping rate κn must satisfy Awκw=Anκn, with Ax the respective cross section areas of both tubes.

As PDs are very narrow, we must take into account that only part of the cross section surface inside the PD is available to a particle of size α. Additionally, a subtle problem lies in the determination of Rw, as it is impossible to create a space filling packing with cylinders. To solve both issues, we rewrite Equation 16 to explicitly contain cross section surfaces. We then replace An with A~~n to accommodate hindrance effects and we replace Aw by 1/ρ. We also ajust PD length: l~=l+2α and L=L-2α. At the same time, we adjust f(σ) to match a triangular distribution of the simplified PDs by using A=1.62 and B=1.36 (Berezhkovskii et al., 2006), which produces the hexagonal cytoplasmic column shape. This yields:

(19) τ=12D[L~2+l~2+2D(l~κn+L~κw)+l~L~(A~~nρ+1A~~nρ)].

We similarly adjust κw:

(20) κw=4ρDHc(α/Rn)Rnf(ρA~n),

where Hc(λ) is the hindrance factor for cylindrical pores (see Materials and methods). In the same fashion, we also adjust κn.

We then invert the relation τ=L22D+L2Peff, where we write Peff for the effective wall permeability (Makhnovskii et al., 2009), to obtain Peff=L2τ-L2/D. With this, we can compute fih=Peff/(ρΠ(α)), where Π(α) is calculated using the same PD geometry. To validate the choice of boundary placement underlying the calculations above, we also calculated the MFPT over two PD passages, that is by shifting the reflecting boundary to the middle of one cell further. This resulted in a 4-fold increase of τ and L2 and hence in exactly the same Peff.

To assess whether the desmotubule has a large impact on fih, we further adapt Equation 19 by replacing A~~n by our desmotubule corrected A~~n, except in f(σ). Additionally, we multiply f(σ) by ξ=(R~n2-R~dt2)/R~n2. Numerical calculations in a simple trapping setup confirm the validity of reducing f(σ) proportional to the area occupied by the desmotubule whilst calculating σ based on the outer radius alone (Appendix 4—figure 1 and Appendix 4). This is in agreement with results for diffusion towards clusters of traps in 3D (Makhnovskii et al., 2000). By the same reasoning, we introduced a hindrance factor in κw. Finally, we adjust the hindrance factors to a slit geometry as before. This results in:

(21) τ=12D[L~2+l~2+L~/ρ+l~A~~n2RnH(2α/(Rn-Rdt))ξf(ρA~n)+l~L~(A~~nρ+1A~~nρ)].

To investigate the effect of different PD distributions, we used all relevant pairs of A and B in f(σ) for different regular trap distributions as given in Berezhkovskii et al. (2006). As Aw is calculated implicitly from 1/ρ, no other adjustments were necessary.

Correction factor fih for pit fields

Request a detailed protocol

For computing fih in pit fields, we used a two step approach similar to computing fih including DT as described above. A similar approach is also followed for the sub-nano channel model. In this calculation, a single pit field is modelled as a number of PDs on a triangular (or square) grid with a centre-to-centre distance d between nearest neighbours. We then calculate the pit radius, Rpit as the radius of the circle that fits the outer edges of the PD entrances. In the trivial case of one PD per ‘pit’, Rpit=Rn. For larger numbers of PDs per pit, see Table 1. For this calculation, individual PDs are modelled as straight cylindrical PDs with radius Rn. We calculate a τ based on circular traps with radius Rpit and a reduced efficiency based on the fraction of the pit that is occupied by the circular PDs. We accordingly adjust κw,pit and τ,pit:

(22) κw,pit=4ρDHc(α/Rpit)Rpitξf(ρπR~pit2),

where p is the number of PDs per pit and ξ=pR~n2/R~pit2 is the fraction of available pit area that is occupied by available PD area, and

(23) τ=12D[L~2+l~2+L~/ρ+l~pA~~n2RpitH(α/Rpit)ξf(ρA~pit)+l~L~(pA~~nρ+1pA~~nρ)].

In these equations, ρ is the total PD density. In our graphs, we either keep ρ constant while increasing p to investigate the effect of clustering, resulting in a pit density ρpits of ρ/p, or keep ρpits constant to investigate the effect of (repeated) PD twinning. As a default, we used d = 120 nm based on distances measured from pictures in Faulkner et al. (2008) of basal cell walls of Nicotinia tabacum leaf trichomes. To verify our calculations, we compared them with a single step calculation with large circles only, that is with radius Rpit and density ρ/p. As results in 3D suggest that for strongly absorbing clusters, the outer radius and cluster density dominate the diffusion (survival time) process (Makhnovskii et al., 2000), this should produce a lower bound to fih. In terms of PDs, this regime applies if a particle that reaches a pit field also has a high probability of entering in it. Indeed, the values calculated with the two step method above were similar and somewhat larger than with the simple large patch method, showing that our computation method is reasonable.

Only a relatively small fraction of the pits is occupied by the PD entrances (5–10% when modelled as circles with Rn = 14 nm and 3–7% with Rn = 12 nm.). Consequently, this approach may become inaccurate when Rpit gets too large. We indeed found instances where fih,pits was larger than fih,singlePDs. In those cases, Rpit was in the order of dpit/4 or larger. We assume that in those cases, the clusters are so close, that the clustering has only minor impact on fih, and fih is better estimated by the calculation for single PDs.

Computing required densities or α¯ with default model

Request a detailed protocol

Numbers in Table 2 are computed based on forward computation of P(α) given ρ, α¯, corresponding Rn and other parameters with increments of 0.1 PD/μm2 (ρ) or 0.01 nm (α¯ etc.) and linear interpolation between the two values that closest match the target P(α). This yields an error of less than 0.0001 μm/s on P(α). We use α = 0.5 nm for CF. The method for computing P(CF) using the unobstructed sleeve (default) model is described throughout the main text. PDinsight, the python program used for computing all values in Table 2, Appendix 5—table 1, Figure 8B and Table 1 is available as supporting material.

Appendix 1

List of mathematical symbols

DefaultDescription
α0.5 nm for CFParticle size
α¯Maximum particle size that fits through a (model) PD
AcCytosolic sleeve cross section area in the central region
A~cCytosolic sleeve cross section area in the central region adjusted for particle size (=steric hindrance only)
A~~cCytosolic sleeve cross section area in the central region after cross sectional averaging (=full hindrance; Ac~~=H(λ)Ac)
AnCytosolic sleeve cross section area in the neck region
CConcentration
ΔCConcentration difference (over the PD channel)
D162 μm2/s for CFParticle size dependent diffusion constant
d11 nm3/sDiffusion constant for particle of unit radius ×1 nm (dummy value to illustrate scaling behaviour)
fihCorrection factor for inhomogeneous wall permeability (0fih1)
H(λ)Hindrance factor calculated for a slit pore geometry
Hc(λ)Hindrance factor calculated for a cylindrical pore geometry
λRelative particle size at the respective location (for straight PDs: λ=α/α¯)
ln25 nmNeck length
l~nEffective neck length for a given particle size (l~n=ln+α)
l100, 200, 500 nmTotal PD length
L10 μmCell length
Rc17.5 nmCentral region radius
Rdt8 nmDT radius
R~dtParticle size adjusted DT radius (R~dt=Rdt+α)
Rn12 nmNeck radius
RxCentral region or neck radius, depending on position
RpitPit field radius: the radius of the smallest circle that circumscribes all PDs in the pit field
R~xParticle size adjusted central region/neck radius (R~x=Rx-α)
ρPD density
p1Number of PDs per cluster (‘pit field’)
P(α)Symplasmic permeability (for particles of size α) of the entire cell wall
Π(α)Symplasmic permeability (for particles of size α) of a single PD per unit of cell wall surface, without correction factor fih (Π(α)=P(α)fihρ).
Q(α)Molar flow rate through a single PD (for particles of size α)
QrelMolar flow rate relative to a reference situation (straight PD)
Q~relRescaled Qrel: Q~rel=(Qrel-1)/(l2l~n-1)+1
τrelMean residence time (MRT) inside PD relative to a reference situation (straight PD)
τ~relRescaled τrel: τ~rel=(τrel-1)l2/(2l~n(l-2l~n))
τMean first passage time (MFPT) through a cytoplasmic column of length L with a wall in the middle of it containing one central PD; Used for the calculation of fih

Appendix 2

Estimating systematic error due to homogeneous flux assumption

Molar flow rate

To estimate the error on the calculation of Q(α), we compared our analytical results (Equation 2) to numerical evaluations of the diffusion equations on a 2D cross section of the available surface of the model PD with neck, see Appendix 2—figure 1.

Appendix 2—figure 1
Error of homogeneous flux approximation (all 2D).

(A) C/x from numerical calculations (2D) along a straight line through the middle of the available neck region for different particle sizes. (B) Top: C/x at neck entrance (proportional to the channel flux) from numerical calculations (N; solid red line with asterisks) and from 3-cylinder model with homogeneous flux assumption (T; dashed cyan line with crosses). The 3-cylinder model results in a consistent over estimation of < 4% (bottom). (C-H) Concentration heat maps for available part of the channel, focus on the neck/central region transition. The same color gradient is used for all six graphs. Black isolines are spaced at 1% of the total concentration difference over the channel. Parameters: l = 200 nm, defaults.

MRT

In absence of a DT, we can compute τ,PD analogously to τ. This yields:

(24) τ,PD=12D(4l~n2+(l2l~n)2+2l~n(l2l~n)(A~~cA~~n+A~~nA~~c)+V~c+2V~n2RnH(λn)f(A~n/A~c)),

with V~x the hindrance adjusted volume of central region or single neck region. The expression in brackets is identical to the one in Equation 14, except for the addition of the last term, meaning that τ is an underestimation of the MFPT. Using τ,PD to define an alternative τrel for channels without DT, τ,rel, suggests that τrel is an underestimate of at least approximately 7–9% for Rc= 17.5 nm and l = 100 nm (Appendix 2—figure 2). This factor saturates between 1.35 and 1.40 for all relevant particle sizes in the limit of unrealistically large Rc. For larger l, the factor has the same maximum values, but these are approached slower. Hence, the error is less for realistic Rc (e.g. 4–6% for l = 200 nm and 2–3% for l = 500 nm.)

Appendix 2—figure 2
Comparison of τrel as used in the main text against a τ,PD based calculation (τ,rel).

l = 100 nm, Rn = 12nm, Rdt = 8 nm, ln = 25 nm.

Appendix 3

Scaling behaviour of Qrel and τrel

Rescaling of Qrel and τrel is a way to collapse our understanding of the processes into simpler curves. The local flux is inversely proportional to the available cross section, motivating the ratio A~~c/A~~n of the dilated channel as a rescaling factor for the x-axis. Using the limit for Qrel it is possible to almost completely collapse the curves for different particle sizes for a single l,ln combination (Appendix 3—figure 1B). For rescaling of τrel, we use its behaviour for large Rc:

(25)τrel(α,Rc)A~~c(l2(l~n+α))(2(l~n+α))l2A~~n,(Rc)(26)2R~c2l~n(l2l~n)l2(R~n2R~dt2),(Rc).
Appendix 3—figure 1
Scaling of Qrel and τrel through the PD.

(A) Maximum possible increase of Q. Line type indicates particle size (solid: α0, dashed: α= 0.5 nm, sparse dashed: α = 1 nm, dash-dotted: α = 1.5 nm). (B) Curves for different α almost collapse with when plotted as function of A~c/A~n and rescaled from the minimal value of 1 to Qrel,max. The curves for different l do not collapse (blue: l = 100 nm, red: l = 200 nm, cyan: l = 500 nm). (C) Curves for τrel collapse (for different l and α) with rescaling function τ~re=(τrel-1)l2/(2l~n(l-2l~n)). Vertical lines in B,C correspond to A~~c/A~~n for different particle sizes (line types as in A). Parameters: ln = 25 nm, Rn = 12 nm, Rdt = 8 nm.

For large Rc, τrel becomes proportional to R~c2. From this we derived a rescaling factor for τrel with f~c the fraction of PD length occupied by the central region (adapted for particle size), that collapses the curves for τrel for all α and l (Appendix 3—figure 1C). The rescaling factor for the x-axis, A~~c/A~~n, increases faster for larger particles. The reason is that A~~n decreases relatively faster with particle size than A~~c, which becomes intuitively clear from Figure 2C,D. This difference explains why the curves for the largest particles are on top prior to rescaling (Figure 4). The τrel-rescaling factor implies that the MRT increases fastest if the central region occupies approximately half of the length of the channel. With our choice of a constant ln = 25 nm, this occurs at a wall thickness of l 100 nm.

Appendix 4

Numerical calculations for trapping rate of annular trap

Seen from the cytoplasmic bulk, the PD orifice has the shape of an annulus (‘ring’). For our calculations above, however, we only have published trapping rates for circular traps. We tested two options for selecting an equivalent circular trap and trapping rate. For this, we numerically solved diffusion equations in a box with a single trap in the middle of one face and a source opposite to it (Appendix 4—figure 1). In the x and y direction, we used periodic boundary conditions reflecting a periodic array of traps. In the z direction we fixed the concentration on side of the domain (‘source plane’) and used a radiating boundary condition with a mixed rate κ(x,y)=k(x,y)D on the other side (‘target plane’). We chose the trapping rate proportional to the diffusion constant, as the flux and molar flow rate through single PDs (Q(α)) are proportional to D (Equation 2). For PDs, the target plane contained the ‘front view’ of a single channel: an annulus with inner radius Rdt, outer radius Rn and surface area Aann. Within this annulus k(x,y) was set to unity (k(x,y)|PD = 1/μm) and 0 outside. For the corresponding homogeneous target plane κhom=AannAtotal/μm. At the grid level, the pixels at a boundary of the annulus had κ proportional to the fraction of their surface falling inside the annulus. The reference flux was computed analytically, exploiting that within each plane at a given distance z from the target plane, the concentration is the same everywhere. This allows for a trivial mapping to a 1D system. These 3D numerical calculations were performed using the Douglas method for 3D alternating direction implicit diffusion (Douglas and Gunn, 1964). To save computation time we used the analytically calculated reference profile as an initial condition for all calculations. We found that the annular patches gave the same result as circles with radius Rn and k(x,y)|circ=AannAcirc/μm, that is that the outer radius of the patch was most important (Appendix 4—figure 1), In line with results for diffusion towards clusters of traps in 3D (Makhnovskii et al., 2000). These calculations support the choice of trapping rate in the calculation of fih including DT.

Appendix 4—figure 1
Correction factor for annular trap shape.

Comparison of discrete and homogeneous permeability. (A) Setup. A homogeneous source is located at a distance h from a trapping plane with either a single trap or homogeneous absorbency. The periodic boundary conditions in the other directions make that the traps are effectively spaced on a square grid with distance d between traps. (B) Relative average net flux in direction of trap (component along the box) for PDs with DT (‘ann’: red) and circular channels with the same apparent surface (‘sCirc’: blue) and same outer radius, but trapping rate decreased according to surface ratios (‘circ’: cyan), h = 300 nm. For both, the flux is compared with a homogeneous trap with rate κ=AannAtotalD/μm. Error bars (nearly invisible) represent the quality of the numerical estimate by minimum and maximum possible values.

Appendix 5

Computing required densities or α¯ with sub-nano model

Required densities or α¯ for the sub-nano model are computed with PDinsight similar to the default model. Additionally, we use the following approximations. For computing ρsn based on a given α¯, we compute a density multiplier based on nc (Equation 7) and Figure 5. This implicitly assumes that fih is not affected by the different PD density. For example, for α¯= 2.5 nm, a 16.8 (Equation 7) × 1.63 (Figure 5) / 9 (channels per PD; Olesen, 1979; Liesche and Schulz, 2013) = 3.04 times as high density would be sufficient. These numbers are the same for straight PDs of any length. For calculating ρsn,neck, we assume that the sub-nano channel structure only occurs in the neck region, with an unobstructed sleeve with Rc=Rdt+2α¯ in the central part. Using a homogeneous flux assumption around the transition between both parts, the factors x reduce to x(2l~+(1-2l~)/x)/l. Similarly, ρgate is computed by assuming a 1 nm thick sub-nano channel structure at both ends of the PD.

As fih is affected by Rn and Rn values get quite large in our calculations with ρ fixed, we follow a different approach for computing α¯sn and Rn,sn based on a given ρ. We use forward calculations based on nine cylindrical channels in a PD, with the trapping rate κw adjusted with an outer radius Rn that would fit all nine channels surrounding the DT.

(27) κw9=4ρDHc(α/Rn)Rnξf(ρπR~n2),

where A~n is the surface area per cylindrical channel and ξ=9(α¯-α)2R~n2 is the fraction of the enveloping circle that is occupied by the nine channels. For sufficiently small α¯, the nine circular channels and minimal protein spacers (at least 1 nm wide) all fit while touching the DT. In that case: Rn=Rn=Rdt+2α¯. With Rdt= 8 nm, this is possible up to α¯ 3.4 nm. For larger α¯, the spacer requirement determines the outer radius of the composite of 9 channels and Rn=α¯+α¯+s/2sin(π/9), where s=1 nm is the spacer width. Rdt does not occur in this equation, because the cylindrical channels can no longer (all) touch the DT.

Appendix 5—table 1
Parameter requirements for reproducing measured P(CF) values based on the default unobstructed sleeve model (also in Table 2) and the sub-nano channel model.

A: Required density (ρ) given α¯ and corresponding neck radius (Rn). B: Required α¯ and corresponding Rn given ρ. For P (CF) = 25μm/s, also values for a 2x, 3x and 4x increase of ρ are computed. This is done both for a uniform increase of the density (p=1) and for (repeated) twinning (p>1) from a uniform starting density (indicated in bold). A, B: The + sign at Rn indicates that the stated Rn is too narrow to fit nine sub-nano channels that touch a DT with Rdt= 8 nm. Models used for calculating required densities: ρ, Rn: default unobstruced sleeve model, ρsn, Rn,sn: sub-nano channel model, ρsn,neck: sub-nano channel model restricted to the neck regions (ln= 25 nm), and ρgate: 1 nm thick structures at both PD entrances locally similar to the sub-nano channel model.

A
P(CF) (μm/s)α¯ (nm) Rn (nm)ρ (PD/μm2)ρsnρsn,neckρgateρsn,neckρgate
62.01233.5139.987.836.7
2.51322.769.146.424.1
3.01416.840.929.117.5
3.414.813.829.121.614.2
3.515+13.227.020.213.6
4.016+10.719.215.010.9
8.52.01247.2197.1123.751.7
2.51332.097.465.434.0
3.01423.757.741.024.7
3.414.819.441.130.520.1
3.515+18.538.128.519.1
4.016+15.027.021.115.4
B
P(CF) (μm/s)ρ α¯Rn (nm)α¯sn*Rn,sn(nm)no spacersRn,sn 1 nmspacers
6104.216.35.218.4+21.8
133.515.14.617.3+19.6
8.5105.218.46.020.1+25.2
134.416.85.418.7+22.5
1101.511.02.613.213.2
131.310.62.412.7 12.7
C
P(CF) (μm/s)ρpα¯Rn (nm)α¯sn*Rn,sn(nm)no spacersRn,sn 1 nmspacers
25 10 1 10.5 28.9 10.0 28.0+ 40.7
2016.621.27.122.2+29.4
27.222.57.523.1+31.1
3015.118.15.919.8+24.6
35.619.26.320.6+26.3
4014.216.45.218.4+21.8
44.617.25.519.0+23.1
13 1 8.8 25.6 8.8 25.5+ 35.9
2615.619.16.320.6+26.2
26.020.06.621.2+27.3
3914.316.65.218.5+22.1
34.617.25.519.0+23.1
5213.615.14.617.3+19.6
43.815.54.817.6+20.4
  1. *: α¯sn is calculated using Rn that allows for 1 nm spacers. C: p is the number of PDs per pit. This table was generated using PDinsight.

Appendix 6

PDinsight

PDinsight is written in python 3. If available, it uses the numpy module, but does not strictly depend upon that. The program has different modes for computing the parameter requirements for a given P(α) and related quantities. The different modes and the relevant parameters are controlled from a parameter file (default: parameters.txt). A graphical user interface (GUI) is provided to help the user create parameter files and run PDinsight. The GUI is written using TKinter, which is included in standard installation of python.

For electron microscopists, who typically have access to many ultrastructural parameters, but often do not know P(α), the mode computeVals will be useful. This computes the expected P(α) values when taking all parameters at face value. Comparison with the sub-nano channel model is possible. In principle, all model parameters must be defined, but missing parameters may be explored using lists of possible values or left at a default (e.g., cell length L and a triangular distribution of PDs), as these have little influence on P(α). If estimates of PD density and distribution are missing, the mode computeUnitVals, which computes Π(α), could be useful. In this mode, comparison with the sub-nano channel is impossible.

In tissue level experiments, P(α) is typically measured, but not all ultrastructural parameters will be known. It is likely that PD density ρ and radius (Rn, assuming straight channels) are poorly known. In this case, mode computeRnDensityGraph will be useful, or a combination of modes computeDens and computeAperture and a number of guesses for Rn/α¯ or ρ, respectively. Additional poorly known parameters can be explored as suggested above. If uncertain, it is strongly advised to explore PD length l (≈ cell wall thickness). For thick cell walls (l 200 nm), it may be worth exploring the effects of increased central radius (Rc>Rn). This is currently only possible in modes computeVals and computeUnitVals.

Major modes

The core of PDinsight is computing effective permeabilities (P(α)) for symplasmic transport based on all model parameters mentioned in the manuscript. This is in mode computeVals. The same computations are used in other modes, which compute the requirements for obtaining a given target value (or set of values) of P(α). In mode computeDens, required densities (ρ) are computed for given values of maximum particle size α¯ (and other parameters). In mode computeAperture, required apertures, given as α¯ as well as neck radius Rn, are computed for given values of ρ (and other parameters). In mode computeRnDensityGraph, Rn,ρ curves are computed that together yield a target P(α). The corresponding values of α¯ are also reported. These curves can be visualized using any plotting program. Mode sensitivityAnalysis computes so-called elasticities (normalized partial derivatives) around a given set (or sets) of parameters. These elasticities tell how sensitive calculated values of P(α),Π(α) and constituents like fih are on the parameters involved.

Auxillary modes

In computing P(α), correction factor fih is automatically included. For specific cases such as modelling studies, however, it may be useful to calculate fih separately. For this purpose, several modes exist for exploring inhomogeneity factor fih: computeFih_subNano (function of α¯; also output values for sub-nano model), computeFih_pitField_dens (function of ρ), computeFih_pitField_xMax (function of α¯) and computeTwinning (function of ρpits, cluster density).

By default, computations are performed for the unobstructed sleeve model. Most computations can also be performed for the sub-nano channel model (Liesche and Schulz, 2013; Comtet et al., 2017). Using switch compSubNano, values for the sub-nano model are also computed.

Graphical user interface

: The GUI to PDinsight is written to facilitate the creation of parameter files and also has a button to run PDinsight directly based on the parameters displayed. In contrast to the generic parameter file, the GUI only shows fields for parameters that are actually used in a specific run mode. The first step of using the gui is to select the run mode (choose from: computeVals, computeUnitVals, computeDens, computeAperture, computeRnDensityGraph and sensitivityAnalysis), followed by ‘load default parameters’. This overwrites any user input from previous modes. All required parameters are shown as text entry fields, with radio buttons for the relevant options. Basic validation occurs on the fly (valid input type, etc). When done, click ‘Run PDinsight’ to write a parameter file and run the program or ‘Export parameter file’ to write a parameter file only. Note that switching modes requires clicking ‘load default parameters’, that is a fresh start.

Additional information can be obtained by clicking the ‘Info’ button and, for certain parameters, clicking on the parameter name. If additional information is available, the mouse cursor changes into a question mark.

Command line usage

(linux): python PDinsight_vXX.py PARAMETERFILE.

(windows): first open a command prompt window (e.g., by searching for ‘cmd’ in the search bar) and go to the directory where PDinsight is located. Then type: py PDinsight_vXX.py PARAMETERFILE

Using the GUI from the command line (linux): python PDinsightGUI_vYY.py or (windows): py PDinsightGUI_vYY.py The correct version of PDinsight_vXX.py should also be available in the current directory.

In the above, XX and YY should be replaced by the respective version numbers.

By default, all outputs are tab separated files (.tsv). The other option is comma separated (.csv), which can be set through the parameter outputType.

The latest version of PDinsight + documentation can be downloaded from https://github.com/eedeinum/PDinsight. DOI: 10.5281/zenodo.3536704

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
    Identification of a developmental transition in plasmodesmatal function during embryogenesis in Arabidopsis thaliana
    1. I Kim
    2. FD Hempel
    3. K Sha
    4. J Pfluger
    5. PC Zambryski
    (2002)
    Development 129:1261–1272.
  39. 39
  40. 40
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
  81. 81
  82. 82
  83. 83

Decision letter

  1. Christian S Hardtke
    Senior Editor; University of Lausanne, Switzerland
  2. Dominique C Bergmann
    Reviewing Editor; Stanford University, United States
  3. Johannes Liesche
    Reviewer; Northwest A&F University, China
  4. Michael Knoblauch
    Reviewer; Washington State University, United States
  5. Valentin Couvreur
    Reviewer; Université catholique de Louvain, Belgium

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

Plasmodesmata (PD) are highly regulated channel-like structures between plant cells that enable movement of diverse molecules including metabolites, proteins and viruses. Long standing questions include how the PD is structured and how this structure enables selectivity that changes during plant growth and upon stress.

In this manuscript, the authors present a mechanistic model to investigate the relation between the geometry of PD, molecular diffusion rate across PD, and molecule residence time in the PD. Key achievements are the scaling of single PD models to whole cell-cell interfaces, providing valuable insight on the design principles of PD and their distributions, as well as providing a tool to predict flux based on structural parameters. This is especially important for internal plant tissues that are not accessible for many experimental approaches. The major novelty of the model is that it is based on a realistic description of PD geometry and mostly relies on measurable parameters. It is a valuable companion to the direct measurements of symplasmic fluxes using molecular reporters and new microscopy technologies allow the fast acquisition of datasets about PD geometry, but that cannot be exploited to their full potential without modelling tools. The release of an open source version of their model constitutes a substantial service to the plant biology community.

Decision letter after peer review:

Thank you for submitting your article "From plasmodesma geometry to effective symplastic permeability through biophysical modelling" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Christian Hardtke as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Johannes Liesche (Reviewer #1); Michael Knoblauch (Reviewer #2); Valentin Couvreur (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Through discussions among the reviewers it was clarified that further experiments would improve the manuscript, but were not strictly necessary for the work to be a strong advance for the field, so long as a number of textual changes were made give a clear context for this work in light of previous studies. For this reason, I am compiling the extensive notes from the three reviewers as a list below. Please consider these previous data when writing a revision and a response.

Overall, it will be essential to demonstrate that the model has broad applicability, but to recognize that it was not yet tested in many situations.

1) Since the model was only applied to a single interface, and reviewer one seems to think the data basis seems to be handled loosely, please see comments from reviewer 1 to address this issue.

2) A more careful evaluation of the current literature, including the discussion of previous models by various authors, is needed.

3) The authors should consider if their model could be even more ambitious, for example by including evaluation of neck length and non-straight shapes.

Reviewer 1:

Introduction paragraph four: The paper by Liesche and Schulz, 2013, is not presented here accurately. This paper compares three different models of PD anatomy. One of them is a cytosolic sleeve model that shares quite some similarity with the model proposed here. The sub-nano channel model was only developed for a very specific question: whether PD can be constricted enough to enable the filter effect that has been ascribed to them as part of the polymer-trap mechanism for phloem loading. According to the hypothesis, sucrose should be able to diffuse through the PD at the bundle-sheath-to-phloem interface, but the slightly larger raffinose should not. It was the aim of that paper to test what kind of PD architecture could be compatible with this mechanism.

Introduction paragraph five: There is no principle limit to using fluorophores of different sizes. Terry and Robards, 1987, actually used tracers of all sizes. Only because fluorescein, especially as carboxyfluorescein diacetate (CFDA) it is widely used, because it crosses membranes, allowing for non-invasive observations.

Modern tracer-based approaches using confocal microscopy can yield very similar results. Compare Rutschow et al., 2011 and Liesche and Schulz, 2013. See also Liesche and Schulz (2012, Plant Physiology), who compared permeabilities across plant species and different cell-cell interfaces.

Photoactivation and photobleaching approaches are not time consuming. Quite the opposite. Application of tracer, sample preparation and imaging can easily done within half an hour.

Introduction paragraph six: I suggest formulating with greater care. Photoactivation and photobleaching approaches as well as GFP transport studies have been very valuable tools for assessing PD function.

The incompatibility here could have very different reasons. A better citation would be Liesche et al., 2019, which actually tries to compare structure-based modeling with functional data. As mentioned above, the Liesche and Schulz, 2013, paper had a different objective.

General comment to the Introduction: I advise a more nuanced view at the current state of the field and approaches carried out in the past. Also, I recommend using the Liesche et al., 2019 paper in Plant Physiology as starting point because it actually is the first systematic comparison of PD structure and function. Moreover, this publication clearly demonstrates the need for better models.

Figure 1: The effective symplasmic permeability in a file of cells also depends on cell size and properties of the cytosol.

Final paragraph of the Introduction: Python is a start, but how about making a standalone app? Or an add-in/macro for Excel? Something that any molecular biologist (and especially students) can use! Last year I got quite a bit of experimental data on phloem transport and I wanted to check it against theoretical models that have been described as easy to use. My basic Matlab skills were not enough to make it work. In my opinion, if you are serious about getting a lot of scientists to try the model, you have to provide a truly easy-to-use solution.

Results subsection “Outline of the model”: Many of the equations used in this model were applied to calculate PD transport for the first time by Liesche and Schulz, 2013. Maybe that should be acknowledged somewhere in this paragraph or the Introduction.

Figure 3: What is the reasoning behind plot B? How can I see from the plot that D always has a big influence (which, of course, is not unexpected)? All values of permeability are very low.

Results paragraph five: density jump is not a term often used in PD literature. You could consider using concentration potential.

Subsection ““Necked” PDs increase molecular flux in thicker cell walls”: I find this heading slightly misleading. Shouldn't it be "dilated PDs increase molecular flux"? Because Rn is the same in both cases, meaning that the neck doesn't change. Or am I just misled by the drawing in Figure 4?

First paragraph of subsection ““Necked” PDs increase molecular flux in thicker cell walls”: Maybe write "electron microscopy" instead of "evidence".

Figure 4: Again, "dilated" seems more appropriate than "necked". Panel C has a different y axis scale from all the other panels. The legend in panel C is confusing as it is missing an α. I recommend either legends in all panels or description in the figure caption.

Figure 4 and accompanying text: Is Rn = 12nm in the 'necked' PD and 17.5 in case of narrow PD? This should be formulated clearly, both in the figure and the text. The drawing above panel A gives the impression that they have the same Rn.

Paragraph three of the same section: In this case the length of the neck region is still 1/3rd of the total PD length, right? It might be worthwhile to also check the effect of variation in the length of the necked region. Ideally, this should also be one parameter of the python tool, because it can be estimated from TEM images.

Subsection “The desmotubule increases PD transport and changes the dependence on particle size”: I also find this heading misleading. My first intuition would be to compare a PD with DT to a PD without DT.

In the same subsection what exactly does the assumption on selection towards size-selectivity mean? That a plant would make the PD as narrow as possible, depending on the molecule that is supposed to pass? This assumption might be valid for relatively big molecules, but probably not for small ones (below 1-2nm). Except for the very special PDs in the polymer trapping species mentioned earlier. I recommend discussing this assumption in more detail.

Where exactly does the "see methods" refer to? I don't find a section on DT-related calculations in the methods.

Paragraph two of subsection “Clustering of PDs in pit 1elds reduces effective symplastic permeability”: Wasn't Rn = 12 nm in Figure 4?

In paragraph three of the same section: Rpit is not in the list of mathematical symbols. Moreover, it is not clear why the radius of a pit field should be comparable to the radius of individual PD, since the number and size of PD within a pit field can vary.

In the final paragraph of the same section: How should this be measurable? As a correlation of clustering with wall thickness?

Figure 7: In panel A, is the distance between PD or the radius of the pit field constant?

Paragraph three of subsection “Application of the model to compute effective permeability for fluorescein derivatives”: With a density of 5 (as specified by Zhu et al., see comments below), it seems like the value for Rn should be 22 nm.

Paragraph four: Why is a density of 10 μm2 used here instead of the value from Zhu et al. (5.4 μm2)?

Also didn't a previous paragraph show that central dilation ("necking") is especially effective for long PDs? How does that fit with the results presented here?

Table 1: PD density seems to start at 10 μm2. This value is higher than the 5.4 μm2 specified by Zhu et al. (Table 2). Their value should be considered a max value as the area that was analyzed by Rutschow et al. is reaching into the mature root zone already. Mature zone has much lower density values (0.62 μm2, Table 3 in Zhu et al.).

In paragraph five the model seems incompatible with the low-concentration H2O2 data. At a density of 5.4 μm2, Rn would need to increase to over 40 nm according to Figure 8B. This seems unrealistic. An increase in PD number can also not be expected because of the short treatment period of 2h.

“We also compared the results obtained with our model and the sub-nano channel model reported before (Liesche and Schulz, 2013).”: This makes very little sense as this model was put forward only for a very specific case, as mentioned above. It would be better to compare with the cytosolic slit model that was also described in that paper.

Discussion paragraph two: Again, this model should not be compared here or it should be clearly stated that it was not developed to be applied in this way. Instead the authors should compare their model to the cytosolic sleeve model of Liesche and Schulz (2013) or the 'pure diffusion' model by Doelger et al., 2014 or the model for 'diffusion through simple PD' included in Ross-Elliott et al., 2017.

“For example, sucrose moves symplastically from bundle sheet cells (BSC) to intermediary cells (IC), where it is polymerized into the larger oligomers raffinose and stachyose, that do not diffuse back in detectable amounts”: It should be mentioned that this only happens in certain species (active symplasmic phloem loaders) with the Cucurbits as the most prominent example.

Discussion paragraph four: The Liesche et al., 2019, paper does not analyze correlation of PD length with permeability at the BSC-phloem interface of active symplasmic loaders. It analyzes correlations of PD anatomical parameters (including length) with permeability across species. I don't understand why a correlation of PD length with permeability would be expected, as the study generally finds permeability to depend to a much higher degree on PD diameter. Please also note that it is very unlikely that a bulk flow was overlooked in that study as the measurements were performed on detached leaves. Indeed, this has been tested by Liesche and Schulz, 2012 who compared full photosynthesising leaf and cut-out tissues and found no difference in PD permeability in tobacco.

Discussion paragraph six: please include 'at the root unloading zone'. Funnel PDs are not found at other phloem interfaces.

“Applying our model for diffusion as a sole driver of symplastic transport can indeed explain experimentally observed measurements of effective symplastic permeability for CF, but only with somewhat wider PDs/neck regions or several fold higher PD densities than usually measured by EM.”: This statement seems to contradict itself. In addition, I highly recommend application of the model to additional interfaces. Liesche et al., 2019, for example, found a model that matches observations of permeability for the bundle-sheath-phloem interface, but not for the bundle sheath-mesophyll interface. In the present case, it should be noted that PD density and dimensions were determined for sand-grown plants (Zhu et al., 1998), whereas effective diffusivity was measured for plate-grown plants (Rutschow et al., 2011). I don't know how growth conditions affect PDs, but it has been previously shown that there are big differences, for example in gene expression, between soil- and plate-grown plants. I mention this just to illustrate the need for additional tests.

“Our model can also explain the effect on permeability after treatment with high and low concentrations of H2O2 in Rutschow et al., 2011.”: Please reconsider this statement regarding low H2O2 in light of comments above.

I am very skeptical about the conclusion. Stress perception, signaling, and two rounds of twinning, all within 2h seem unrealistic. I also didn't find evidence for such a rapid multiplication of PD numbers in the two citations provided here. Moreover, 30 PD μm2 would be an extremely high number, that, to my knowledge, has never been observed at a "regular" cell-cell interface.

Discussion paragraph seven, “Despite these deviations, comparing our model to the sub-nano channel model, we found that the latter requires roughly twice as high PD densities to produce the same permeability values P (CF)” and following sentences: Again, this is not an appropriate comparison. Please compare to other models listed in a comment above.

https://doi.org/10.7554/eLife.49000.sa1

Author response

Through discussions among the reviewers it was clarified that further experiments would improve the manuscript, but were not strictly necessary for the work to be a strong advance for the field, so long as a number of textual changes were made give a clear context for this work in light of previous studies. For this reason, I am compiling the extensive notes from the three reviewers as a list below. Please consider these previous data when writing a revision and a response.

Overall, it will be essential to demonstrate that the model has broad applicability, but to recognize that it was not yet tested in many situations.

We now explicitly discuss the need for additional testing and its expected benefits in the manuscript. To illustrate broad applicability, we have analyzed three additional experiments where pure diffusive transport was expected (an epidermis only experiment in Rutshow et al., 2011, from which we already analyzed the other experiments; additional data obtained from Johannes Liesche (reviewer 1), intended for a follow up manuscript analyzing even more data. See responses to reviewer 1 and Table 1 for more details. As argued in the introduction, pure diffusive symplasmic transport is likely to be the relevant form of symplasmic transport in many developmentally relevant processes/conditions. Additionally, we have added a short section to the Discussion describing what would be needed to include advection/flow while conserving the unique features of our model. Finally, to allow more researchers to join in this broader application effort, we have put considerable effort in building a graphical user interface to the python program, which, according to testers, was makes it much easier to use.

1) Since the model was only applied to a single interface, and reviewer one seems to think the data basis seems to be handled loosely, please see comments from reviewer 1 to address this issue.

We assume that the suggestion that “the data basis seems to be handled loosely” is mostly based on the comments by reviewer 1. As described there, we disagree on the remark that we should have used a 2-fold lower PD density. The density of ~5PD μm2 applies to the epidermis. About 2 times higher densities are reported for more interior layers, which comprise the bulk of the tissue in the experiment we analyze. For consistency, we now also analyzed the second measurement from the Rutschow paper, which is limited to the epidermis and where lower permeabilities are reported from the experiment. In our analysis, the different permeabilities reported in Rutschow et al., 2011 and the different densities for different interfaces reported in Zhu et al., 1999 are all consistent.

To test for broader applicability, we have also applied the model (with slight addition of functionality) to two data sets originally acquired in the context of Liesche et al., 2019: mesophyll-mesophyll interfaces in poplar and Cucurbita max. As stated in the general response, we have included the results of this analysis in a separate document for internal use in the reviewing process, but prefer not to include them in the current manuscript but make them part of a larger analysis / collaborative follow up project, 1) to avoid a potential conflict of interests and 2) because their inclusion would not cancel the suggestion that further testing remains highly valuable (and is an important reason for building a community tool).

2) A more careful evaluation of the current literature, including the discussion of previous models by various authors, is needed.

We have incorporated a discussion of previous models by various authors in the Introduction. See responses to reviewer 1 and general comment on the Introduction for details. We also spend more words in the discussion comparing our model to various other models using a similar PD architecture.

3) The authors should consider if their model could be even more ambitious, for example by including evaluation of neck length and non-straight shapes.

We have included a figure on the impact of varying neck length between 15 – 50 nm (Figure 4—figure supplement 2). Variation of neck length (and different neck lengths/diameters on either side of the PD) has now been included in PDinsight in the process of analyzing the additional data. The suggestion of non-straight shapes is certainly interesting and envisioned as a future extension of PDinsight, but careful implementation of this feature is not feasible within the two months given for resubmission and, moreover, would add yet another topic to an already dense manuscript.

Reviewer 1:

Introduction paragraph four: The paper by Liesche and Schulz, 2013, is not presented here accurately. This paper compares three different models of PD anatomy. One of them is a cytosolic sleeve model that shares quite some similarity with the model proposed here. The sub-nano channel model was only developed for a very specific question: whether PD can be constricted enough to enable the filter effect that has been ascribed to them as part of the polymer-trap mechanism for phloem loading. According to the hypothesis, sucrose should be able to diffuse through the PD at the bundle-sheath-to-phloem interface, but the slightly larger raffinose should not. It was the aim of that paper to test what kind of PD architecture could be compatible with this mechanism.

We originally cited Liesche and Schulz, 2013 for this model architecture as the first to model this architecture. To avoid confusion, we have emphasized the specific context used in Liesche and Schultz, 2013, which is the same context as in Comtet et al., 2017, now also cited for using the sub-nano channel architecture.

The sub-nano architecture is now part of a section that describes the generally used architectures, i.e., a cytosolic slit without further localized obstructions and the sub-nano channel structure. Liesche and Schultz, 2013 is cited in both categories.

Introduction paragraph five: There is no principle limit to using fluorophores of different sizes. Terry and Robards, 1987, actually used tracers of all sizes. Only because fluorescein, especially as carboxyfluorescein diacetate (CFDA) it is widely used, because it crosses membranes, allowing for non-invasive observations.

We agree that the text was insufficiently clear at this point. We have added a few words emphasizing that the size limitation relates to the non-invasive techniques and additionally emphasizing the wider size range of microinjection/bombardment methods.

Modern tracer-based approaches using confocal microscopy can yield very similar results. Compare Rutschow et al., 2011 and Liesche and Schulz, 2013. See also Liesche and Schulz (2012, Plant Physiology), who compared permeabilities across plant species and different cell-cell interfaces.

We assume the Liesche and Schulz, 2013, is the following 2012 manuscript: https://doi.org/10.1111/j.1365-2818.2011.03584.x We have added a reference to this manuscript stating values in the same range (and some lower ones).

Photoactivation and photobleaching approaches are not time consuming. Quite the opposite. Application of tracer, sample preparation and imaging can easily done within half an hour.

Introduction paragraph six: I suggest formulating with greater care. Photoactivation and photobleaching approaches as well as GFP transport studies have been very valuable tools for assessing PD function.

We have rephrased the sentence to make it correct: time consuming now only relates to the generation of transgenic lines with cell specific fluorescent probes (as originally intended).

The incompatibility here could have very different reasons. A better citation would be Liesche et al., 2019, which actually tries to compare structure-based modeling with functional data. As mentioned above, the Liesche and Schulz, 2013, paper had a different objective.

Agreed. Reference changed.

General comment to the Introduction: I advise a more nuanced view at the current state of the field and approaches carried out in the past. Also, I recommend using the Liesche et al., 2019 paper in Plant Physiology as starting point because it actually is the first systematic comparison of PD structure and function. Moreover, this publication clearly demonstrates the need for better models.

In response to this general comment, we have made the following changes.

1) Including a more extensive overview of modelling efforts, including manuscripts focussing on flow through PDs. In this, we have classified the models based on the general PD architecture used: an unobstructed sleeve surrounding a desmotubule (the majority) and rarer geometries:

– Funnel shape (only Ross-Elliott et al., 2017) in the phloem unloading zone.

– The sub-nano channel architecture (Liesche and Schulz, 2013 and Comtet et al., 2017), where we also mention the specific context of size selectivity for small molecules.

2) Rephrasing the limitations of modern techniques.

3) Towards the end of the Introduction, we have used Liesche et al., 2019 to illustrate the need for better models, as suggested above.

In processing all reviewer comments, we also realized another issue that hampers the coupling of tissue level and ultrastructural findings: the tissue level measurements are not always reported in (SI) units that can be compared to computed effective symplasmic permeabilities. We have now mentioned this in the Discussion.

Figure 1: The effective symplasmic permeability in a file of cells also depends on cell size and properties of the cytosol.

As we use it, effective symplasmic permeability is a property of the cell wall only. We use cell length L (which is proportional to cell size in a file of rectangular cells) to compute the correction factor fih, but this completely cancels out of the equations. Consequently, P(α) is fully determined by PDs/wall properties (and particle size). Therefore, cell size and cytosol properties are not included in Figure 1. To emphasize the above, we have now explicitly stated this aim when introducing the model setup.

Final paragraph of the Introduction: Python is a start, but how about making a standalone app? Or an add-in/macro for Excel? Something that any molecular biologist (and especially students) can use! Last year I got quite a bit of experimental data on phloem transport and I wanted to check it against theoretical models that have been described as easy to use. My basic Matlab skills were not enough to make it work. In my opinion, if you are serious about getting a lot of scientists to try the model, you have to provide a truly easy-to-use solution.

To increase the ease of use, we have developed a graphical user interface (GUI) that aids the user in creating a correct parameter file (all required numbers etc can be entered in the GUI, which also performs basic verification) and then allows the user to run PDinsight by clicking a button. The GUI is written in TKinter, which is included in any standard Python installation. We have also explored the option of an excel macro/add-in, but we concluded this would be too inflexible and hard to maintain. In the future, we will continue to think about ways to improve the program, working both on easier user interaction and adding/adjusting core functionality to the demands of real data. We already implemented additional features to optimally analyze the two Liesche et al., 2019 derived data sets mentioned above.

Compared to the version uploaded for review, we have made the following improvements:

– GUI for easier access.

– Improved documentation (currently available as appendix to the manuscript).

– Default parameter files for each of the modes, containing only those parameters that actually are used in the specific computation.

– Added a mode for computing a ``unit permeability’’ (\Pi(\α) = “permeation constant of a single PD” = equation 4, in the manuscript), useful for comparisons if density information is lacking.

– Removed the dependency on numpy (if installed, numpy is still used), so the user does not have to install any python modules to make the program work.

Results subsection “Outline of the model”: Many of the equations used in this model were applied to calculate PD transport for the first time by Liesche and Schulz, 2013,. Maybe that should be acknowledged somewhere in this paragraph or the Introduction.

At a conceptual level, this is now acknowledged in the Introduction when mentioning the use of an unobstructed sleeve architecture.

Many equations follow from the geometrical assumptions, but the use of hindrance factors is a point where models differ a lot. So, for this crucial point we added a reference to Liesche and Schulz, 2013 when introducing the use of hindrance factors.

Figure 3: What is the reasoning behind plot B? How can I see from the plot that D always has a big influence (which, of course, is not unexpected)? All values of permeability are very low.

The reasoning behind plot B is to show that the particle’s diffusion constant D has a very large impact on permeability and explains a large part of its size dependence. This is visible in plot B by seeing similar slopes and, for sufficiently large Rn, also values among the curves. For the largest Rn, there is only a 3-fold difference left between 0.1 nm and 2 nm particles (in B), whereas the difference is 50-60-fold in the more realistic plot A. We have added a line in the text describing Figure 3B with this explanation: “For example, at $R_n=R_c$, the 50+ fold difference between $\α=0.1$ nm and $\α=2$ nm is reduced to a 3-fold difference.”

Re low values: All values are in arbitrary units. Therefore, their absolute value is meaningless. We have added a sentence to emphasize this.

Results paragraph five: density jump is not a term often used in PD literature. You could consider using concentration potential.

We have opted for “concentration difference”.

Subsection ““Necked” PDs increase molecular flux in thicker cell walls”: I find this heading slightly misleading. Shouldn't it be "dilated PDs increase molecular flux"? Because Rn is the same in both cases, meaning that the neck doesn't change. Or am I just misled by the drawing in Figure 4?

Agreed, changed to “A dilated central region increases molecular flux in thicker cell walls”.

First paragraph of subsection ““Necked” PDs increase molecular flux in thicker cell walls”: Maybe write "electron microscopy" instead of "evidence".

Agreed, changed as suggested.

Figure 4: Again, "dilated" seems more appropriate than "necked". Panel C has a different y axis scale from all the other panels. The legend in panel C is confusing as it is missing an α. I recommend either legends in all panels or description in the figure caption.

Figure 4 and accompanying text: Is Rn = 12nm in the 'necked' PD and 17.5 in case of narrow PD? This should be formulated clearly, both in the figure and the text. The drawing above panel A gives the impression that they have the same Rn.

The neck diameter in all cases is 12 nm as the drawing on top suggests. We have clarified this in the text and by adding a sentence to the beginning of the caption of Figure 4.

We have changed the top to “dilated” and added α to the legend (also in Figure 4—figure supplement 1). We have not changed the y axes of all panels, because that would reduce visibility. This change also affects the parameter Qrel and \taurel.

Paragraph three of the same section: In this case the length of the neck region is still 1/3rd of the total PD length, right? It might be worthwhile to also check the effect of variation in the length of the necked region. Ideally, this should also be one parameter of the python tool, because it can be estimated from TEM images.

Added a figure (Figure 4—figure supplement 2) exploring the effects of neck region length. Neck length can already be entered as a parameter in PDinsight.

Subsection “The desmotubule increases PD transport and changes the dependence on particle size”: I also find this heading misleading. My first intuition would be to compare a PD with DT to a PD without DT.

Changed heading, emphasizing that PDs with a given maximum particle size (~SEL) are compared. Note that we briefly investigate the effect of DT removal in Figure 8CD in the context of the Rutschow et al., 2011 experiments.

In the same subsection what exactly does the assumption on selection towards size-selectivity mean? That a plant would make the PD as narrow as possible, depending on the molecule that is supposed to pass? This assumption might be valid for relatively big molecules, but probably not for small ones (below 1-2nm). Except for the very special PDs in the polymer trapping species mentioned earlier. I recommend discussing this assumption in more detail.

Size selectivity is meant in the context of the model: control over the maximum particle size that can pass. This could yield a requirement that at a certain developmental stage, only particles of, say, <2.5 nm radius should be allowed to pass, which in our model translates directly to a 5 nm difference between Rdt and Rn, or, for cylindrical channels without DT, a radius of 2.5 nm.

To avoid confusion (the same point was also raised by reviewer 2), we have rephrased the sentence and removed the reference to natural selection.

Where exactly does the "see methods" refer to? I don't find a section on DT-related calculations in the methods.

Reference removed. The full calculation is included in the main text.

Paragraph two of subsection “Clustering of PDs in pit 1elds reduces effective symplastic permeability”: Wasn't Rn = 12 nm in Figure 4?

Corrected.

In paragraph three of the same section: Rpit is not in the list of mathematical symbols. Moreover, it is not clear why the radius of a pit field should be comparable to the radius of individual PD, since the number and size of PD within a pit field can vary.

Added Rpit to the list. Also added a few words to emphasise that Rpit is much larger than Rn (except for the trivial case with a single PD in a “pit field” that no biologist would call pit field).

In the final paragraph of the same section: How should this be measurable? As a correlation of clustering with wall thickness?

The aim of the sentence was to emphasize that for realistic Rn, the effects of Figure 6B,D would be too small to be experimentally measurable in any way, followed by the observation that PD length and density have more impact on fih when PDs are clustered into pit fields. We have rephrased the sentence to avoid this confusion.

We agree that measuring fih experimentally is currently impossible, as the discrepancies in permeability values from bleaching/photoactivation experiments and calculated from EM observations are too large with current techniques. This may remain so for a long time, and even if technically possible, the question is who would be interested in measuring this correction factor independently.

Figure 7: In panel A, is the distance between PD or the radius of the pit field constant?

Distance between PDs. Added a few words to clarify. The cartoons are indeed suggesting otherwise, which is hard to prevent without using cartoons that are too large.

Paragraph three of subsection “Application of the model to compute effective permeability for fluorescein derivatives”: With a density of 5 (as specified by Zhu et al., see comments below), it seems like the value for Rn should be 22 nm.

Required Rn would be 16.5, 16.9, or 18.3 nm, depending on the assumed PD density (5.42/ 5 / 4 PD μm2, respectively), a target value of P=3.3 um/s and otherwise the same parameters as in our main text. We have now also included the calculations for this second experiment in the main text and Table 1. Initially, we had not included this second experiment in our manuscript, as for geometrical reasons the interpretation of the single cell bleach is more complicated than the tissue level bleach. This is illustrated by the fact that Rutschow et al. chose to neglect any coupling in the radial and tangential directions (detailed in the supplementary text of Rutschow et al).

Paragraph four: Why is a density of 10 μm2 used here instead of the value from Zhu et al. (5.4 μm2)?

Zhu et al. Table 2 reports 5.4 density values for the epidermis whereas the reported values for tangential walls of inner cortex (12.58), outer cortex (9.05) and (and similar to vascular (9.92)) were in the range of densities used here (10 μm2 -13 μm2). The quoted targets of P=6-8.5 um/s were based on the tissue level bleaches, which involve multiple cell layers and hence include only a limited fraction of epidermal cells (Figure 1 Rutschow et al).

Rutschow et al. also report measurements for single epidermal cells, where the value of 5.42 would be applicable, possibly weighted with the radial value of 3.33 μm2. For this measurement (Figure 3 Rutschow et al), they report P=3.3 +/- 0.8 um/s. (also see previous comment).

To clarify this issue we have explicitly included the exact measurements of Zhu et al.. We have also included calculations for the epidermal experiment, there using a density of 5 (or 5.42) μm2. The predictions are remarkably consistent.

Also didn't a previous paragraph show that central dilation ("necking") is especially effective for long PDs? How does that fit with the results presented here?

Unfortunate phrasing, rephrased. The point is that central dilation can increase permeability (most effectively in long PDs), but only to a limited extent (a straight channel with radius Rc). The Rutschow values cannot be reproduced by central dilation.

Table 1: PD density seems to start at 10 μm2. This value is higher than the 5.4 μm2 specified by Zhu et al. (Table 2). Their value should be considered a max value as the area that was analyzed by Rutschow et al. is reaching into the mature root zone already. Mature zone has much lower density values (0.62 μm2, Table 3 in Zhu et al.).

See above: the value of 5.4 μm2 is not applicable to the experiment we analyzed. As we have now also analyzed the epidermal experiment (with a target effective symplasmic permeability of 3.3 um/s), we have also added 5 and 5.42 μm2 to Table 1B.

The Figure 1 experiment in Rutschow et al. was performed in a 50 μm zone at ~200 μm from the QC. Cell shapes seem not elongated yet (eg, see Figure 3 in the Rutschow paper), so likely before the “transition zone” (after Sabatini).

The measurements for immature tissue in Zhu et al. were performed at 100 / 150 μm from the QC, depending on cell type / interface (not 100% clear from the manuscript which was used for each measurement). Their measurements for mature tissue were performed at 1200 – 1250 μm from the QC.

Taken together, while we appreciate the concern, it seems likely that the densities in the Rutschow experiment are closer to the immature than the “mature” densities in Zhu et al. We have, therefore, not added the low “mature” densities to Table 1.

Another way to look at the data is that the root zone used in the Rutschow experiment has to have “immature” PD densities if PD transport is predominantly diffusive in this zone/direction, because otherwise absurd PD densities and/or Rn would be required for reproducing the measured permeabilities assuming only diffusive transport.

In paragraph five the model seems incompatible with the low-concentration H2O2 data. At a density of 5.4 μm2, Rn would need to increase to over 40 nm according to Figure 8B. This seems unrealistic. An increase in PD number can also not be expected because of the short treatment period of 2h.

See comments above: a density of ~10 μm2 is more likely and the speed of the sink-source transition in (tobacco) leaves and corresponding changes in PD structure indicate that a density increase could occur in a short window of time, although such a massive increase in 2 hours has not been reported under natural circumstances.

That said, explaining these data seem to require a dramatic change at the PD level of whatever kind. We have altered the text and more explicitly in the Discussion, to emphasize this point.

“We also compared the results obtained with our model and the sub-nano channel model reported before (Liesche and Schulz, 2013).”: This makes very little sense as this model was put forward only for a very specific case, as mentioned above. It would be better to compare with the cytosolic slit model that was also described in that paper.

We agree that comparison to other cytosolic sleeve models is important, so we have added this in the Discussion. The comparison with the sub-nano channel architecture is valuable because it is also used by others (Comtet et al., 2017), and as a (pessimistic) lower bound for the performance of a sleeve with spokes. We have rephrased the text and removed the reference to Liesche and Schulz at this point to emphasize that the comparison here is about different possible geometries, not specific instances of the model. Also see the next comment.

Discussion paragraph two: Again, this model should not be compared here or it should be clearly stated that it was not developed to be applied in this way. Instead the authors should compare their model to the cytosolic sleeve model of Liesche and Schulz, 2013, or the 'pure diffusion' model by Doelger et al., 2014 or the model for 'diffusion through simple PD' included in Ross-Elliott et al., 2017.

We have kept a short comparison between sub-nano channel and unobstructed sleeve _architecture_, which we subsequently use to state that the sub-nano calculations provide a (pessimistic) estimate for the effect of spokes in the cytosolic sleeve. To avoid the suggestion that the sub-nano channel model was developed to explain the Rutschow et al. data, we compare “architectures” and do not reference Liesche et al., 2013 here.

In the two paragraphs above, we have included a comparison with other models: the three straight unobstructed sleeve models mentioned by the reviewer 629-632) and the central dilation geometry used by Blake, 1978 (who only studies flow, not diffusion).

Differences in model PD geometry are also detailed in the Introduction.

“For example, sucrose moves symplastically from bundle sheet cells (BSC) to intermediary cells (IC), where it is polymerized into the larger oligomers raffinose and stachyose, that do not diffuse back in detectable amounts”: It should be mentioned that this only happens in certain species (active symplasmic phloem loaders) with the Cucurbits as the most prominent example.

Mention added.

Discussion paragraph four: The Liesche et al., 2019, paper does not analyze correlation of PD length with permeability at the BSC-phloem interface of active symplasmic loaders. It analyzes correlations of PD anatomical parameters (including length) with permeability across species. I don't understand why a correlation of PD length with permeability would be expected, as the study generally finds permeability to depend to a much higher degree on PD diameter. Please also note that it is very unlikely that a bulk flow was overlooked in that study as the measurements were performed on detached leaves. Indeed, this has been tested by Liesche and Schulz, 2012 (Plant Physiology) who compared full photosynthesising leaf and cut-out tissues and found no difference in PD permeability in tobacco.

Considering the issues mentioned above as well as other issues that came up during subsequent discussion among us authors, we have decided to remove the sentence and reference here.

Discussion paragraph six: please include 'at the root unloading zone'. Funnel PDs are not found at other phloem interfaces.

Added.

“Applying our model for diffusion as a sole driver of symplastic transport can indeed explain experimentally observed measurements of effective symplastic permeability for CF, but only with somewhat wider PDs/neck regions or several fold higher PD densities than usually measured by EM.”: This statement seems to contradict itself. In addition, I highly recommend application of the model to additional interfaces. Liesche et al., 2019, for example, found a model that matches observations of permeability for the bundle-sheath-phloem interface, but not for the bundle sheath-mesophyll interface. In the present case, it should be noted that PD density and dimensions were determined for sand-grown plants (Zhu et al., 1998), whereas effective diffusivity was measured for plate-grown plants (Rutschow et al., 2011). I don't know how growth conditions affect PDs, but it has been previously shown that there are big differences, for example in gene expression, between soil- and plate-grown plants. I mention this just to illustrate the need for additional tests.

We agree that testing more interfaces is important. We have added a section emphasizing the need for and benefits of additional testing following the discussion of this experiment. Moreover, we have modified the discussion structure, in the process of merging both points where the Rutschow et al., 2011 experiments are discussed. As mentioned, we now have analyzed both Rutschow experiments, one involving an entire root section, the other only the epidermis. Additionally, we have analyzed data from two mesophyll-mesophyll interfaces in poplar and Cucurbita max that were originally acquired in the context of Liesche et al., 2019, but not published in the detailed form that we needed. Results of this analysis are included in a separate document for the reviewers/editors, but are not meant to be included in the current manuscript to avoid a possible conflict of interest.

“Our model can also explain the effect on permeability after treatment with high and low concentrations of H2O2 in Rutschow et al., 2011.”: Please reconsider this statement regarding low H2O2 in light of comments above.

We have rephrased the statement such that high H2O2 can easily be explained, whereas low H2O2 is much harder. We have also rewritten the remainder of the section and omitted the references mentioned below.

I am very skeptical about the conclusion. Stress perception, signaling, and two rounds of twinning, all within 2h seem unrealistic. I also didn't find evidence for such a rapid multiplication of PD numbers in the two citations provided here. Moreover, 30 PD μm2 would be an extremely high number, that, to my knowledge, has never been observed at a "regular" cell-cell interface.

We agree that the 30+ density is indeed very high. Then again, regular cell-cell interfaces are not regularly exposed to 2 hours of H2O2 signaling. The only way to find out is to prepare similarly treated tissue for EM investigation. If the density is indeed dramatically increased, this should be relatively easy to see. We have added this suggested experiment to the text.

Additionally, we have added some more skepticism in the discussion of the low H2O2 results. The conclusion of (limited) evidence in those papers is based on the following observations: Roberts et al., 2001: formation of complex PDs coincides with the leaf sink-source transition in tobacco leaves. Particularly the walls of pallisade and spongy parenchyma show sharp transitions in the occurrence of branched PDs (Figure 5 of Roberts et al.). At a macroscopic level, the sink-source transition progressed over the leave with speeds of 100 um/h (small leaves, avg 4.5 cm), 310 um/h (medium leaves, avg 7 cm), 600 um/h (large leaves, avg 10 cm) (measured using AtSuc2-GFP), i.e., in all cases involving many cells per hour. → With a speed of multiple cell walls per hour, PD architecture might be changed from simple to predominantly complex, which involves a large increase in the number of PD entrances at the transformed PDs.

Fitzgibbon et al., 2013: “With young detached leaves (leaves 1 and 2), we noticed that the conversion of simple to complex PDs occurred rapidly. New complex PDs could be visualized within 1 h of leaf detachment and were formed continuously at an average rate of 0.18 PDs per cell per hour (Figure 2G). In the areas observed, the overall numbers of PDs increased two- to threefold within 48 h.” → formation of complex PDs (from simple), which may have more than 2 channels, is possible within 1 h. The H2O2 treatment would have to synchronize the process and induce it in parallel at many PDs in the same interface, as the “natural” rate of 0.18 PDs per cell per hour is far insufficient.

We thought that this consideration was too detailed and too long for something that produces a mere hint that the required density might be possible (but certainly no strong evidence) to be included in the text. Therefore, we have deleted the references to the above papers.

Discussion paragraph seven, “Despite these deviations, comparing our model to the sub-nano channel model, we found that the latter requires roughly twice as high PD densities to produce the same permeability values P (CF)” and following sentences: Again, this is not an appropriate comparison. Please compare to other models listed in a comment above.

We have restructured the Discussion, in the process reducing emphasis on the difference between sub-nano and unobstructed sleeve architecture, as well as including more comparison with other unobstructed sleeve models including the ones above. As mentioned before, we think that the comparison with the sub-nano channel architecture is valuable, so not removed.

https://doi.org/10.7554/eLife.49000.sa2

Article and author information

Author details

  1. Eva E Deinum

    Mathematical and statistical methods (Biometris), Wageningen University, Wageningen, Netherlands
    Contribution
    Conceptualization, Software, Formal analysis, Funding acquisition, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    eva.deinum@wur.nl
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-8564-200X
  2. Bela M Mulder

    1. Living Matter Department, Institute AMOLF, Amsterdam, Netherlands
    2. Laboratory of Cell Biology, Wageningen University, Wageningen, Netherlands
    Contribution
    Supervision, Methodology, Writing—original draft, Writing—review and editing
    Contributed equally with
    Yoselin Benitez-Alfonso
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8620-5749
  3. Yoselin Benitez-Alfonso

    Centre for Plant Science, University of Leeds, Leeds, United Kingdom
    Contribution
    Conceptualization, Supervision, Writing—original draft, Writing—review and editing
    Contributed equally with
    Bela M Mulder
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9779-0413

Funding

European Molecular Biology Organization (ASTF 105 - 2012)

  • Eva E Deinum

Nederlandse Organisatie voor Wetenschappelijk Onderzoek

  • Bela M Mulder

Engineering and Physical Sciences Research Council (EF/M027740/1)

  • Yoselin Benitez-Alfonso

Leverhulme Trust (RPG-2016-13)

  • Yoselin Benitez-Alfonso

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

The authors thank EM Bayer, W Nicolas, M Glavier and L Brocard from the Bordeaux Imaging Centre, Plant Imaging Platform, ((http://www.bic.u-bordeaux.fr), for providing the electron micrograph included in Figure 1. Also we thank Philip Kirk (Centre for Plant Sciences, Leeds) for critical review of the python program.

Senior Editor

  1. Christian S Hardtke, University of Lausanne, Switzerland

Reviewing Editor

  1. Dominique C Bergmann, Stanford University, United States

Reviewers

  1. Johannes Liesche, Northwest A&F University, China
  2. Michael Knoblauch, Washington State University, United States
  3. Valentin Couvreur, Université catholique de Louvain, Belgium

Publication history

  1. Received: June 3, 2019
  2. Accepted: November 16, 2019
  3. Accepted Manuscript published: November 22, 2019 (version 1)
  4. Accepted Manuscript updated: November 25, 2019 (version 2)
  5. Version of Record published: January 31, 2020 (version 3)

Copyright

© 2019, Deinum et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,239
    Page views
  • 228
    Downloads
  • 2
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cell Biology
    2. Genetics and Genomics
    Weina Zhang et al.
    Tools and Resources
    1. Cell Biology
    2. Developmental Biology
    Xiaofei Bai et al.
    Research Article Updated