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


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.


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.


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.

  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)
P(CF)(μm/s)ρα¯Rn (nm)
P(CF)(μm/s)ρpα¯Rn (nm)
  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.


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:


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)],


(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

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


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:

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.

P(CF) (μm/s)α¯ (nm) Rn (nm)ρ (PD/μm2)ρsnρsn,neckρgateρsn,neckρgate
P(CF) (μm/s)ρ α¯Rn (nm)α¯sn*Rn,sn(nm)no spacersRn,sn 1 nmspacers
131.310.62.412.7 12.7
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
13 1 8.8 25.6 8.8 25.5+ 35.9
  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 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 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 PARAMETERFILE

Using the GUI from the command line (linux): python or (windows): py The correct version of 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 DOI: 10.5281/zenodo.3536704

Data availability

PDinsight can be downloaded from GitHub: Documentation on the use of is included as an appendix to the manuscript with additional information at the head of the example parameter file. More extensive documentation is included with PDinsight on GitHub. PDinsight also has a citable DOI through Zenodo: 10.5281/zenodo.3536704. The PDinsight parameter files used for this manuscript are included as Source code 1.


    1. Kim I
    2. Hempel FD
    3. Sha K
    4. Pfluger J
    5. Zambryski PC
    Identification of a developmental transition in plasmodesmatal function during embryogenesis in Arabidopsis thaliana
    Development 129:1261–1272.

Article and author information

Author details

  1. Eva E Deinum

    Mathematical and statistical methods (Biometris), Wageningen University, Wageningen, Netherlands
    Conceptualization, Software, Formal analysis, Funding acquisition, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing
    For correspondence
    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
    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
    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


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.


The authors thank EM Bayer, W Nicolas, M Glavier and L Brocard from the Bordeaux Imaging Centre, Plant Imaging Platform, ((, 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.

Version 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)


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


  • 2,049
  • 327
  • 27

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

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)

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

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

  1. Eva E Deinum
  2. Bela M Mulder
  3. Yoselin Benitez-Alfonso
From plasmodesma geometry to effective symplasmic permeability through biophysical modelling
eLife 8:e49000.

Share this article

Further reading

    1. Cell Biology
    2. Computational and Systems Biology
    N Suhas Jagannathan, Javier Yu Peng Koh ... Lisa Tucker-Kellogg
    Research Article

    Bats have unique characteristics compared to other mammals, including increased longevity and higher resistance to cancer and infectious disease. While previous studies have analyzed the metabolic requirements for flight, it is still unclear how bat metabolism supports these unique features, and no study has integrated metabolomics, transcriptomics, and proteomics to characterize bat metabolism. In this work, we performed a multi-omics data analysis using a computational model of metabolic fluxes to identify fundamental differences in central metabolism between primary lung fibroblast cell lines from the black flying fox fruit bat (Pteropus alecto) and human. Bat cells showed higher expression levels of Complex I components of electron transport chain (ETC), but, remarkably, a lower rate of oxygen consumption. Computational modeling interpreted these results as indicating that Complex II activity may be low or reversed, similar to an ischemic state. An ischemic-like state of bats was also supported by decreased levels of central metabolites and increased ratios of succinate to fumarate in bat cells. Ischemic states tend to produce reactive oxygen species (ROS), which would be incompatible with the longevity of bats. However, bat cells had higher antioxidant reservoirs (higher total glutathione and higher ratio of NADPH to NADP) despite higher mitochondrial ROS levels. In addition, bat cells were more resistant to glucose deprivation and had increased resistance to ferroptosis, one of the characteristics of which is oxidative stress. Thus, our studies revealed distinct differences in the ETC regulation and metabolic stress responses between human and bat cells.

    1. Cell Biology
    2. Developmental Biology
    Filip Knop, Apolena Zounarova ... Marie Macůrková
    Research Article

    During Caenorhabditis elegans development, multiple cells migrate long distances or extend processes to reach their final position and/or attain proper shape. The Wnt signalling pathway stands out as one of the major coordinators of cell migration or cell outgrowth along the anterior-posterior body axis. The outcome of Wnt signalling is fine-tuned by various mechanisms including endocytosis. In this study, we show that SEL-5, the C. elegans orthologue of mammalian AP2-associated kinase AAK1, acts together with the retromer complex as a positive regulator of EGL-20/Wnt signalling during the migration of QL neuroblast daughter cells. At the same time, SEL-5 in cooperation with the retromer complex is also required during excretory canal cell outgrowth. Importantly, SEL-5 kinase activity is not required for its role in neuronal migration or excretory cell outgrowth, and neither of these processes is dependent on DPY-23/AP2M1 phosphorylation. We further establish that the Wnt proteins CWN-1 and CWN-2 together with the Frizzled receptor CFZ-2 positively regulate excretory cell outgrowth, while LIN-44/Wnt and LIN-17/Frizzled together generate a stop signal inhibiting its extension.