From plasmodesma geometry to effective symplasmic permeability through biophysical modelling
Abstract
Regulation of molecular transport via intercellular channels called plasmodesmata (PDs) is important for both coordinating developmental and environmental responses among neighbouring cells, and isolating (groups of) cells to execute distinct programs. Celltocell 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 multilevel model allows to predict PD dimensions matching measured permeabilities and add a functional interpretation to structural differences observed between PDs in different cell walls.
Introduction
The formation of spatial patterns in plants requires the transport and interaction of molecular signals. This sharing of information coordinates cell fate decisions over multiple cells and the isolation of cell fate determinants within a cell or group of cells on the same developmental path. Small molecules such as sugars, peptides, hormones and RNAs move long and short distances to coordinate cell/organ development (Otero et al., 2016). Celltocell 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 SHORTROOT (SHR), an Arabidopsis thaliana GRAS family transcription factor, that moves from the stele to corticalendodermal 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 (cytoplasmtocytoplasm) 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 beta1,3 glucan polysaccharide synthesized by callose synthases and degraded by β−1,3glucanases (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; BenitezAlfonso et al., 2009; Xu et al., 2012).
Small molecules can move via PD by diffusion (nontargeted 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 (RossElliott 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; BretHarte and Silk, 1994; Jensen et al., 2012; RossElliott 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 (BretHarte and Silk, 1994; Liesche and Schulz, 2013; Dölger et al., 2014; RossElliott 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. (RossElliott 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 socalled ‘subnano 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 9fold 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 DRONPAs (28 kDa), Dendra2 (26 kDa), (photoactive and nonphotoactive) 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 cellspecific 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 multispecies 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 nontargeted symplasmic transport and reveal the potential source of conflicts with ultrastructural measurements. We found that, in this context, our model performed better than the ‘subnano channel model’ (Liesche and Schulz, 2013) referred to above. Our calculations demonstrate that multiscale modelling approaches can integrate results from PD structural dimensions and molecular fluxes and reveal conflicts on these determinations. We, therefore, recommend these should be applied systematically when defining effective symplasmic permeability for a particular tissue/molecule and/or biological context. To facilitate this, we share a python program for computing effective permeabilities from PD geometries as a community resource.
Results
Outline of the model
Our aim is to describe the symplasmic transport properties of a cell wall as an effective wall permeability, that is a single number that could be plugged into tissue/organ level models. For this, we split the transport into two parts: the movement through an individual channel representing a PD and the approach to this channel from the cytoplasmic bulk (Figure 1). This implicitly assumes a homogeneous cytosol. The basic geometrical terminology that we considered in our calculations is introduced in the cartoon PD shown in Figure 1B. An overview of all mathematical symbols is given in Appendix 1.
Obtaining good EM data of PD dimensions is notoriously hard. We therefore opted for a simple geometrical description that allows us to study the effects of PD neck, central region and desmotubule dimensions with as few parameters as possible (see Materials and methods). We modelled a single PD as a 3part 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 ${l}_{n}$ and radius ${R}_{n}$. The central region has radius ${R}_{c}$. Over the whole length, the center of the channel is occupied by a ‘desmotubule’ (DT) modelled as a cylinder of radius ${R}_{dt}$. The part available for diffusive transport, the cytosolic sleeve, is the space between the outer cylinder wall and the DT.
We made the arguably simplest choice of modelling particles as (nonadditive, i.e. not interacting among themselves) hard spheres with radius $\alpha $. 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 $\alpha $ (Figure 2B,C). This socalled 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, $\overline{\alpha}$, 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 $\overline{\alpha}$ to operationalize the SEL concept in a straightforward manner. To avoid confusion, however, we will consistently write $\overline{\alpha}$ when referring to our model.
We introduced rescaled geometrical parameters to account for the reduced available volume in a compact way: ${\stackrel{~}{l}}_{n}={l}_{n}+\alpha $, ${\stackrel{~}{R}}_{c}={R}_{c}\alpha $, ${\stackrel{~}{R}}_{dt}={R}_{dt}+\alpha $ and ${\stackrel{~}{R}}_{n}={R}_{n}\alpha $. With these, the available surface area (Figure 2C) is
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: $\overline{\alpha}=({R}_{n}{R}_{dt})/2$.
Considering pure diffusion without particle turnover inside the PD, particle flux through the channel is described by $\frac{\partial {C}_{xyz}}{\partial t}=D{\nabla}^{2}{C}_{xyz}$, or in steady state: $D{\nabla}^{2}{C}_{xyz}=0$, with ${C}_{xyz}$ 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$, $\nabla {C}_{x}$, is inversely proportional to the available surface area ${A}_{x}$, so $\nabla {C}_{c}={\stackrel{~}{A}}_{n}/{\stackrel{~}{A}}_{c}\nabla {C}_{n}$. The total concentration difference over the PD, $\mathrm{\Delta}C={C}_{l}{C}_{0}$ is accordingly distributed over the channel: $\mathrm{\Delta}C=2\stackrel{~}{{l}_{n}}\nabla {C}_{n}+(l2\stackrel{~}{{l}_{n}})\nabla {C}_{c}$. The steady state molar flow rate $Q(\alpha )$ through each channel is proportional to the entrance cross section: $Q(\alpha )=D{\stackrel{~}{A}}_{n}\nabla {C}_{n}$. Solving these equations for $\nabla {C}_{n}$ leads to:
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 socalled hindrance factors $0\le H(\lambda )<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 slitpore geometry (see Materials and methods). These factors depend on the relative particle size $\lambda $. In our case, $\lambda =2\alpha /({R}_{x}{R}_{dt})$. In the neck region, $\lambda =\alpha /\overline{\alpha}$. For the full expression and behaviour of $H(\lambda )$, see Materials and methods.
As $H(\lambda )$ already includes the effect of steric hindrance between wall and particle, we can adjust Equation 2 by replacing every instance of $\stackrel{~}{{A}_{x}}$ with
For completeness, we note that the simplification of a uniform particle flux along the channel axis is violated near the neckcentral 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, $\mathrm{\Pi}(\alpha )$, through the rule rule steadystate flow rate = permeation constant × concentration difference, yielding
We also defined $\tau $ 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).
Having defined the permeation constant of a single channel, the effective symplasmic permeability of the wall as a whole ($P(\alpha )$, the quantity that can be estimated using tissue level measurements) follows from the definition $J=P\mathrm{\Delta}C$ ($\text{steady state flux}=\text{permeability}\times \text{density jump}$):
with $\rho $, the density of PDs per unit wall area (number/ μm^{2}) and ${f}_{ih}$, a (density dependent) correction factor for the inhomogeneity of the wall ($0<{f}_{ih}<1$). The latter takes into account that the wall is, in fact, only permeable at discrete spots. To calculate ${f}_{ih}$, 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 $\rho \mathrm{\Pi}(\alpha )$ (as described in the Materials and methods).
As expected, $P(\alpha )$ depends on particle size. Two factors underlie this size dependence, which both affect $\mathrm{\Pi}(\alpha )$: hindrance effects, which reduce the space available for particle diffusion, and the fact that the diffusion constant is inversely proportional to particle size: $D={d}_{1}/\alpha $. Figure 3A and (Figure 3—figure supplement 1) show that hindrance effects have the strongest impact for particle sizes close to the maximum $\overline{\alpha}$, whereas the particle diffusion constant always has a large impact Figure 3B. For example, at ${R}_{n}={R}_{c}$, the 50+ fold difference between α = 0.1 nm and α = 2 nm is reduced to a 3fold difference when ignoring the particle size dependence of the diffusion constant.
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 $\overline{\alpha}$). We investigated how both the transport volume (using Equation 2) and transport time ($\tau $ as above) change when the central region is dilated. To compare channels with neck and dilated central region (12 nm $={R}_{n}\le {R}_{c}$) with narrow straight channels (${R}_{n}={R}_{c}=$ 12 nm), we define a relative molar flow rate as ${Q}_{rel}={Q}_{dilated}/{Q}_{narrow}$ and similarly relative ${\tau}_{rel}$ (Figure 4). For a more detailed discussion of ${\tau}_{rel}$ and its computation, see Materials and methods and Appendix 2.
We then investigated how both ${Q}_{rel}$ and ${\tau}_{rel}$ change with increasing central region radius ${R}_{c}$ and how this depends on particle radius $\alpha $ 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 ($\alpha $, dashed lines). In fact, from studying the limiting behaviour of the underlying formulas, we found that ${Q}_{rel}$ is always less than its theoretical maximum $\frac{l}{2{\stackrel{~}{l}}_{n}}$, whereas ${\tau}_{rel}$ ultimately scales quadratically with ${R}_{c}$, and, equivalently, linearly with the surface ratio ${\stackrel{~}{\stackrel{~}{A}}}_{c}/{\stackrel{~}{\stackrel{~}{A}}}_{n}$ (see Appendix 3 and Appendix 3—figure 1). In simpler terms: the benefits of increased transport volume with increasing ${R}_{c}$ saturate, and instead the costs in transport time increases ever faster with further dilation of the central region. This defines a tradeoff between transport volume and transport time with increasing ${R}_{c}$ 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 ${Q}_{rel}$ is possible (Figure 4A–C). Similarly, for any given combination of ${R}_{n}$ and ${R}_{c}$, ${Q}_{rel}$ decreases with increasing ${l}_{n}$ and decreases faster for shorter $l$, whereas ${\tau}_{rel}$ has its maximum at ${\stackrel{~}{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 (${R}_{c}$ = 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 ${R}_{c}$ = 26.4 nm and ${R}_{c}$ = 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 ${R}_{c}$ = 17.5 nm (${R}_{c}^{*}$ in A,B). We found that the molar flow rate ${Q}_{rel}$ increases less than the MRT ${\tau}_{rel}$ when increasing ${R}_{c}$ 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 ${R}_{c}$ = 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 $\overline{\alpha}$. Assuming that control over maximum particle size $\overline{\alpha}$ 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(\alpha )$ is proportional to orifice area ($\approx {A}_{n}$), we first computed ${n}_{c}(\overline{\alpha})$, the number of circular channels that would offer the same ${A}_{n}$ as a single channel with a DT of radius $R}_{dt$ = 8 nm and the same $\overline{\alpha}$:
Figure 5A displays the ${n}_{c}(\overline{\alpha})$ as a function of the maximum particle size. As an example, when $\overline{\alpha}$ = 2 nm, 20 cylindrical channels without DT would be needed to match the orifice surface area of a single channel with DT (with $R}_{dt$ = 8 nm). This number decreases for larger $\overline{\alpha}$. 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 ($\alpha /\overline{\alpha}$). 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).
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’, ${f}_{ih}$ in Equation 6 for the effective symplasmic permeability. Here, we explore how ${f}_{ih}$ depends on all model parameters. To calculate ${f}_{ih}$, 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 ${f}_{ih}$ 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 ${f}_{ih}$ as a function of ${R}_{n}$ and explored its dependence on particle size $\alpha $ (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 $\rho $ (B), wall thickness $l$ (C) and PD distribution in the wall (D). We found that, provided that ${R}_{n}$ is large enough for particles to enter (as indicated by vertical cyan lines in Figure 6—figure supplement 1A), ${f}_{ih}$ is independent of cell length $L$ and particle size $\alpha $ (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 ${f}_{ih}$ also hardly depends on the precise layout of PDs (Figure 6D). Although variations in ${f}_{ih}$ appear larger at low PD densities, for typical ${R}_{n}$ values (for example, 12 nm as in Figure 4) density only has a minor impact (Figure 6B). Finally, we found an increase of ${f}_{ih}$ with increasing PD length $l$, saturating to its theoretical maximum of ${f}_{ih}=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.
Second, we investigated the effect of PDs grouped in small clusters resembling pit fields (see Materials and methods). The average centretocentre 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 ${f}_{ih}$ 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 ${f}_{ih}$ decreases with increasing number of PDs in a pit (and constant total PD density $\rho $). Different from isolated PDs, Figure 7A also reveals that, when grouped in pit fields, there is a strong dependence of ${f}_{ih}$ 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 ${R}_{pit}$ are much larger than the largest ${R}_{n}$ used in Figure 6B. Figure 7B shows that clustering (in this case 7 PDs) increases the dependence of ${f}_{ih}$ 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 ${f}_{ih}$ on PD density. Also the arrangement of PDs in small model clusters affects the degree of dependence ${f}_{ih}$ on $\rho $. In both cases, we observe the steepest dependency of ${f}_{ih}$ on $\rho $ 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).
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(\alpha )$ always increased with the increase in cluster size/PD number (Figure 7D), despite the decrease in ${f}_{ih}$ 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 ${f}_{ih}$ for inhomogeneous wall permeability has only a minor role on $P(\alpha )$. For realistic PD dimensions ($R}_{n$ < 20–25 nm), the additional effect of ${f}_{ih}$ 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, ${f}_{ih}$ is markedly reduced, and PD length and density have a much larger impact on ${f}_{ih}$. 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 nontargeted 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 nonfluorescent 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/μm^{2}. Based on these numbers we assume a PD density of 10–13 PDs/μm^{2} for the tissue level experiment and 5 PDs/μm^{2 }for 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 μm^{2}/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 ($R}_{n$ > 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 ${R}_{dt}$, ${R}_{n}$ 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 $R}_{n$ = 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).
Using the model, we also explored the effect of ‘necked’/‘dilated’ PDs by adding a wider central region to PDs. For a central radius $R}_{c$ = 17.6 nm, the required ${R}_{n}$ to reproduce the tissue level CF permeability values would decrease by perhaps 1 nm or at most 3 nm (for $R}_{c$ = 26.4 nm) considering a PD density in the order of $\rho$ = 10 μm^{–2} (Figure 8C, ${R}_{c}$ 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 H_{2}O_{2} treatment. They found a strong decrease in symplasmic permeability to ≈ 1μm/s after treatment with a ‘high’ H_{2}O_{2} 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 ${R}_{n}$ to 11 nm ($\rho$ = 10 μm^{–2}) or 10.6 nm ($\rho$ = 13 μm^{–2}), resulting in $\overline{\alpha}$ = 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’ H_{2}O_{2} concentration. Reproducing this increase requires a large change at the PD level. At the extremes, an increase of ${R}_{n}$ to approximately 29 nm for $\rho$ = 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 ${R}_{n}$ and $\rho $ would have to increase substantially (Figure 8B). As an extreme hypothesis, we also calculated the effects of complete DT removal. The increases in $P(\alpha )$ that could be obtained this way were by far insufficient to explain the reported effect of mild H_{2}O_{2} 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 subnano channel architecture. Using the subnanochannel architecture, much larger PD densities would be required to achieve the same $P(CF)$: roughly twice as large for $\overline{\alpha}$ = 3.5–4 nm and even larger for smaller $\overline{\alpha}$ (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 H_{2}O_{2} treatment on effective permeability.
Discussion
In this manuscript, we presented a method for calculating effective wall permeabilities for nontargeted, 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’ ${f}_{ih}$ 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 $\overline{\alpha}$ 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 (BenitezAlfonso 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 curvatureinducing 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 (${R}_{c}$), 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 narrowwidenarrow 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 $\overline{\alpha}$, 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 noninvasive techniques, introduces a large uncertainty in SEL measurements and hence $\overline{\alpha}$. Also other biological factors could lead to an underestimation as well as an overestimation of $\overline{\alpha}$. For example, in socalled 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 $\overline{\alpha}$ could occur for nonspherical 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 (BenitezAlfonso et al., 2010) or because membranes and/or cell wall domains around PDs allow for reversible transient modifications in $\overline{\alpha}$ (AbouSaleh et al., 2018). Additionally, molecules could pass in the wake of larger proteins/complexes/structures that modify PDs (e.g., tubuleforming 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 socalled ‘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; RossElliott et al., 2017). All else remaining equal, ‘twinned’ PDs have twice the entrance surface area, which would result in doubling the effective permeability $P(\alpha )$. 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 subchannels (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 ${f}_{ih}$ after twinning suggest that minimizing the distance between subchannels 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 (RossElliott 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. (RossElliott 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 necklike 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 (BenitezAlfonso et al., 2009; BenitezAlfonso et al., 2013; Nicolas et al., 2017b), and similar to ${R}_{c}$ 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 ≈ 2fold 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 sideeffects 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 H_{2}O_{2} in Rutschow et al. (2011). The reduced permeability after high H_{2}O_{2} 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/μm^{2}, a reduction from $\overline{\alpha}$ = 4.2–5.2 nm to $\overline{\alpha}$ = 1.5 nm would be required, see Table 2B). The strongly increased permeability after low H_{2}O_{2} treatment, however, is harder to explain. With a single parameter change, the model predicts either a very wide PD aperture of $\overline{\alpha}$ = 8.8–10.5 nm, or a ±4fold 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 ${R}_{n}$ because plants could respond to preEM 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 ${R}_{n}$ (AbouSaleh et al., 2018; Yan et al., 2019; Amsbury and BenitezAlfonso, 2019). For a passive transport mechanism, the elastic energy required for these reversible deformations would have to be in the order of a few ${k}_{B}T$ 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 ${f}_{ih}$ 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 (BretHarte and Silk, 1994; Liesche and Schulz, 2013; Dölger et al., 2014; RossElliott et al., 2017). Using typical PD dimensions and no clustering, inhomogeneity factor ${f}_{ih}$ 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 sin^{2} 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(\alpha )$ (Appendix 2), but results in parameters that are directly measurable on EM images.
Comparing the unobstructed sleeve architecture to the subnano 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 subnano 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 ${f}_{ih}$ 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 electrontomography (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 protocolSimilar 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 protocolHindrance factors $H(\lambda )$ 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.
This choice is supported by the steric hindrance prefactor that is included in $H(\lambda )$ (Dechadilok and Deen, 2006). This $\mathrm{\Phi}(\lambda )=1\lambda $ is the same as the ratio of available to full surface area ${\stackrel{~}{A}}_{x}(\alpha )/{A}_{x}$. For cylindrical channels, that is reference channels in Figure 5 and the regular PDs after DT removal, we use
for $\lambda <0.95$ and the asymptotic approximation by Mavrovouniotis and Brenner (1988),
otherwise, as suggested by Dechadilok and Deen (2006).
Relative molar flow rate and MRT
Request a detailed protocolFor assessing the impact of the neck constriction on PD transport, we defined two relative quantities: ${Q}_{rel}={Q}_{dilated}/{Q}_{narrow}$ and ${\tau}_{rel}={\tau}_{dilated}/{\tau}_{narrow}$ (Figure 4, Appendix 3—figure 1). Using Equation 2 for $Q(\alpha )$, ${Q}_{rel}$ is well defined:
For ${\tau}_{rel}$ we first needed an expression for $\tau $ itself. Ideally, this would be a MFPT, which could calculated in a way similar to ${\tau}_{\parallel}$ in the calculation of ${f}_{ih}$, using a narrowwidenarrow setup. These calculations, however, critically depend on trapping rates at the narrowwide 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 ${f}_{ih}$, 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(\alpha )$ and defined $\tau $ 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 ${R}_{c}$. In our definition of ${\tau}_{rel}$, the concentration difference cancels from the equation, solving the problem:
This method of computing ${\tau}_{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 protocolTo compute ${f}_{ih}$, 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) ${\tau}_{\parallel}$ through a simplified PD and a column of cytoplasm surrounding it. We then convert ${\tau}_{\parallel}$ to an effective wall permeability and compare the result with the uncorrected effected permeability computed using Equation 6 for the simplified PD geometry and ${f}_{ih}=1$.
As a simplified PD, we use a narrow cylindrical channel of length $l$ and radius ${R}_{n}$, 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 ${R}_{w}$ with a hexagonal column with cross section area ${A}_{w}=1/\rho $ 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:
where
is a trapping rate to map the 3D setup onto a 1D diffusion problem. In this,
is a function that monotonically increases from 0 to infinity as $\sigma $, the fraction of the wall occupied by the circular PDs, increases from 0 to 1. $f(\sigma )$ 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 ${\kappa}_{n}$ must satisfy ${A}_{w}{\kappa}_{w}={A}_{n}{\kappa}_{n}$, with ${A}_{x}$ 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 $\alpha $. Additionally, a subtle problem lies in the determination of ${R}_{w}$, 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 ${A}_{n}$ with ${\stackrel{~}{\stackrel{~}{A}}}_{n}$ to accommodate hindrance effects and we replace ${A}_{w}$ by $1/\rho $. We also ajust PD length: $\stackrel{~}{l}=l+2\alpha $ and $L=L2\alpha $. At the same time, we adjust $f(\sigma )$ 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:
We similarly adjust ${\kappa}_{w}$:
where ${H}_{c}(\lambda )$ is the hindrance factor for cylindrical pores (see Materials and methods). In the same fashion, we also adjust ${\kappa}_{n}$.
We then invert the relation ${\tau}_{\parallel}=\frac{{L}^{2}}{2D}+\frac{L}{2{P}_{eff}}$, where we write ${P}_{eff}$ for the effective wall permeability (Makhnovskii et al., 2009), to obtain ${P}_{eff}=\frac{L}{2{\tau}_{\parallel}{L}^{2}/D}$. With this, we can compute ${f}_{ih}={P}_{eff}/(\rho \mathrm{\Pi}(\alpha ))$, where $\mathrm{\Pi}(\alpha )$ 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 4fold increase of ${\tau}_{\parallel}$ and ${L}^{2}$ and hence in exactly the same ${P}_{eff}$.
To assess whether the desmotubule has a large impact on ${f}_{ih}$, we further adapt Equation 19 by replacing ${\stackrel{~}{\stackrel{~}{A}}}_{n}$ by our desmotubule corrected ${\stackrel{~}{\stackrel{~}{A}}}_{n}$, except in $f(\sigma )$. Additionally, we multiply $f(\sigma )$ by $\xi =({\stackrel{~}{R}}_{n}^{2}{\stackrel{~}{R}}_{dt}^{2})/{\stackrel{~}{R}}_{n}^{2}$. Numerical calculations in a simple trapping setup confirm the validity of reducing $f(\sigma )$ proportional to the area occupied by the desmotubule whilst calculating $\sigma $ 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 ${\kappa}_{w}$. Finally, we adjust the hindrance factors to a slit geometry as before. This results in:
To investigate the effect of different PD distributions, we used all relevant pairs of $A$ and $B$ in $f(\sigma )$ for different regular trap distributions as given in Berezhkovskii et al. (2006). As ${A}_{w}$ is calculated implicitly from $1/\rho $, no other adjustments were necessary.
Correction factor ${f}_{ih}$ for pit fields
Request a detailed protocolFor computing ${f}_{ih}$ in pit fields, we used a two step approach similar to computing ${f}_{ih}$ including DT as described above. A similar approach is also followed for the subnano channel model. In this calculation, a single pit field is modelled as a number of PDs on a triangular (or square) grid with a centretocentre distance $d$ between nearest neighbours. We then calculate the pit radius, ${R}_{pit}$ as the radius of the circle that fits the outer edges of the PD entrances. In the trivial case of one PD per ‘pit’, ${R}_{pit}={R}_{n}$. For larger numbers of PDs per pit, see Table 1. For this calculation, individual PDs are modelled as straight cylindrical PDs with radius ${R}_{n}$. We calculate a ${\tau}_{\parallel}$ based on circular traps with radius ${R}_{pit}$ and a reduced efficiency based on the fraction of the pit that is occupied by the circular PDs. We accordingly adjust ${\kappa}_{w,pit}$ and ${\tau}_{\parallel ,pit}$:
where $p$ is the number of PDs per pit and $\xi =p{\stackrel{~}{R}}_{n}^{2}/{\stackrel{~}{R}}_{pit}^{2}$ is the fraction of available pit area that is occupied by available PD area, and
In these equations, $\rho $ is the total PD density. In our graphs, we either keep $\rho $ constant while increasing $p$ to investigate the effect of clustering, resulting in a pit density ${\rho}_{pits}$ of $\rho /p$, or keep ${\rho}_{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 ${R}_{pit}$ and density $\rho /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 ${f}_{ih}$. 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 $R}_{n$ = 14 nm and 3–7% with $R}_{n$ = 12 nm.). Consequently, this approach may become inaccurate when ${R}_{pit}$ gets too large. We indeed found instances where ${f}_{ih,pits}$ was larger than ${f}_{ih,singlePDs}$. In those cases, ${R}_{pit}$ was in the order of ${d}_{pit}/4$ or larger. We assume that in those cases, the clusters are so close, that the clustering has only minor impact on ${f}_{ih}$, and ${f}_{ih}$ is better estimated by the calculation for single PDs.
Computing required densities or $\overline{\alpha}$ with default model
Request a detailed protocolNumbers in Table 2 are computed based on forward computation of $P(\alpha )$ given $\rho $, $\overline{\alpha}$, corresponding ${R}_{n}$ and other parameters with increments of 0.1 PD/μm^{2} ($\rho $) or 0.01 nm ($\overline{\alpha}$ etc.) and linear interpolation between the two values that closest match the target $P(\alpha )$. This yields an error of less than 0.0001 μm/s on $P(\alpha )$. We use $\alpha$ = 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
Default  Description  

$\alpha $  0.5 nm for CF  Particle size 
$\overline{\alpha}$  Maximum particle size that fits through a (model) PD  
${A}_{c}$  Cytosolic sleeve cross section area in the central region  
${\stackrel{~}{A}}_{c}$  Cytosolic sleeve cross section area in the central region adjusted for particle size (=steric hindrance only)  
${\stackrel{~}{\stackrel{~}{A}}}_{c}$  Cytosolic sleeve cross section area in the central region after cross sectional averaging (=full hindrance; $\stackrel{~}{\stackrel{~}{{A}_{c}}}=H(\lambda ){A}_{c}$)  
${A}_{n}$  Cytosolic sleeve cross section area in the neck region  
$C$  Concentration  
$\mathrm{\Delta}C$  Concentration difference (over the PD channel)  
$D$  162 μm^{2}/s for CF  Particle size dependent diffusion constant 
${d}_{1}$  1 nm^{3}/s  Diffusion constant for particle of unit radius ×1 nm (dummy value to illustrate scaling behaviour) 
${f}_{ih}$  Correction factor for inhomogeneous wall permeability ($0\le {f}_{ih}\le 1$)  
$H(\lambda )$  Hindrance factor calculated for a slit pore geometry  
${H}_{c}(\lambda )$  Hindrance factor calculated for a cylindrical pore geometry  
$\lambda $  Relative particle size at the respective location (for straight PDs: $\lambda =\alpha /\overline{\alpha}$)  
${l}_{n}$  25 nm  Neck length 
${\stackrel{~}{l}}_{n}$  Effective neck length for a given particle size (${\stackrel{~}{l}}_{n}={l}_{n}+\alpha $)  
$l$  100, 200, 500 nm  Total PD length 
$L$  10 μm  Cell length 
${R}_{c}$  17.5 nm  Central region radius 
${R}_{dt}$  8 nm  DT radius 
${\stackrel{~}{R}}_{dt}$  Particle size adjusted DT radius (${\stackrel{~}{R}}_{dt}={R}_{dt}+\alpha $)  
${R}_{n}$  12 nm  Neck radius 
${R}_{x}$  Central region or neck radius, depending on position  
${R}_{pit}$  Pit field radius: the radius of the smallest circle that circumscribes all PDs in the pit field  
${\stackrel{~}{R}}_{x}$  Particle size adjusted central region/neck radius (${\stackrel{~}{R}}_{x}={R}_{x}\alpha $)  
$\rho $  PD density  
$p$  1  Number of PDs per cluster (‘pit field’) 
$P(\alpha )$  Symplasmic permeability (for particles of size $\alpha $) of the entire cell wall  
$\mathrm{\Pi}(\alpha )$  Symplasmic permeability (for particles of size $\alpha $) of a single PD per unit of cell wall surface, without correction factor ${f}_{ih}$ ($\mathrm{\Pi}(\alpha )=\frac{P(\alpha )}{{f}_{ih}\rho}$).  
$Q(\alpha )$  Molar flow rate through a single PD (for particles of size $\alpha $)  
${Q}_{rel}$  Molar flow rate relative to a reference situation (straight PD)  
${\stackrel{~}{Q}}_{rel}$  Rescaled ${Q}_{rel}$: ${\stackrel{~}{Q}}_{rel}=({Q}_{rel}1)/(\frac{l}{2{\stackrel{~}{l}}_{n}}1)+1$  
${\tau}_{rel}$  Mean residence time (MRT) inside PD relative to a reference situation (straight PD)  
${\stackrel{~}{\tau}}_{rel}$  Rescaled ${\tau}_{rel}$: ${\stackrel{~}{\tau}}_{rel}=({\tau}_{rel}1){l}^{2}/(2{\stackrel{~}{l}}_{n}(l2{\stackrel{~}{l}}_{n}))$  
${\tau}_{\parallel}$  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 ${f}_{ih}$ 
Appendix 2
Estimating systematic error due to homogeneous flux assumption
Molar flow rate
To estimate the error on the calculation of $Q(\alpha )$, 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.
MRT
In absence of a DT, we can compute ${\tau}_{\parallel ,PD}$ analogously to ${\tau}_{\parallel}$. This yields:
with ${\stackrel{~}{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 $\tau $ is an underestimation of the MFPT. Using ${\tau}_{\parallel ,PD}$ to define an alternative ${\tau}_{rel}$ for channels without DT, ${\tau}_{\parallel ,rel}$, suggests that ${\tau}_{rel}$ is an underestimate of at least approximately 7–9% for $R}_{c$= 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 ${R}_{c}$. For larger $l$, the factor has the same maximum values, but these are approached slower. Hence, the error is less for realistic ${R}_{c}$ (e.g. 4–6% for $l$ = 200 nm and 2–3% for $l$ = 500 nm.)
Appendix 3
Scaling behaviour of ${Q}_{rel}$ and ${\tau}_{rel}$
Rescaling of ${Q}_{rel}$ and ${\tau}_{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 ${\stackrel{~}{\stackrel{~}{A}}}_{c}/{\stackrel{~}{\stackrel{~}{A}}}_{n}$ of the dilated channel as a rescaling factor for the xaxis. Using the limit for ${Q}_{rel}$ it is possible to almost completely collapse the curves for different particle sizes for a single $l,{l}_{n}$ combination (Appendix 3—figure 1B). For rescaling of ${\tau}_{rel}$, we use its behaviour for large ${R}_{c}$:
For large ${R}_{c}$, ${\tau}_{rel}$ becomes proportional to $\stackrel{~}{R}}_{c}^{2$. From this we derived a rescaling factor for $\tau}_{rel$ with ${\stackrel{~}{f}}_{c}$ the fraction of PD length occupied by the central region (adapted for particle size), that collapses the curves for ${\tau}_{rel}$ for all $\alpha $ and $l$ (Appendix 3—figure 1C). The rescaling factor for the xaxis, ${\stackrel{~}{\stackrel{~}{A}}}_{c}/{\stackrel{~}{\stackrel{~}{A}}}_{n}$, increases faster for larger particles. The reason is that ${\stackrel{~}{\stackrel{~}{A}}}_{n}$ decreases relatively faster with particle size than ${\stackrel{~}{\stackrel{~}{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 ${\tau}_{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 $l}_{n$ = 25 nm, this occurs at a wall thickness of $l\approx$ 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 $\kappa (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(\alpha )$) are proportional to $D$ (Equation 2). For PDs, the target plane contained the ‘front view’ of a single channel: an annulus with inner radius ${R}_{dt}$, outer radius ${R}_{n}$ and surface area ${A}_{ann}$. 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 ${\kappa}_{hom}=\frac{{A}_{ann}}{{A}_{total}}/\mu m$. At the grid level, the pixels at a boundary of the annulus had $\kappa $ 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 ${R}_{n}$ and ${k(x,y)}_{circ}=\frac{{A}_{ann}}{{A}_{circ}}/\mu 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 ${f}_{ih}$ including DT.
Appendix 5
Computing required densities or $\overline{\alpha}$ with subnano model
Required densities or $\overline{\alpha}$ for the subnano model are computed with PDinsight similar to the default model. Additionally, we use the following approximations. For computing ${\rho}_{sn}$ based on a given $\overline{\alpha}$, we compute a density multiplier based on ${n}_{c}$ (Equation 7) and Figure 5. This implicitly assumes that ${f}_{ih}$ is not affected by the different PD density. For example, for $\overline{\alpha}$= 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 ${\rho}_{sn,neck}$, we assume that the subnano channel structure only occurs in the neck region, with an unobstructed sleeve with ${R}_{c}={R}_{dt}+2\overline{\alpha}$ in the central part. Using a homogeneous flux assumption around the transition between both parts, the factors $x$ reduce to $x(2\stackrel{~}{l}+(12\stackrel{~}{l})/x)/l$. Similarly, ${\rho}_{gate}$ is computed by assuming a 1 nm thick subnano channel structure at both ends of the PD.
As ${f}_{ih}$ is affected by ${R}_{n}$ and ${R}_{n}$ values get quite large in our calculations with $\rho $ fixed, we follow a different approach for computing ${\overline{\alpha}}_{sn}$ and ${R}_{n,sn}$ based on a given $\rho $. We use forward calculations based on nine cylindrical channels in a PD, with the trapping rate ${\kappa}_{w}$ adjusted with an outer radius ${R}_{n}^{\prime}$ that would fit all nine channels surrounding the DT.
where ${\stackrel{~}{A}}_{n}$ is the surface area per cylindrical channel and $\xi =\frac{9{(\overline{\alpha}\alpha )}^{2}}{{\stackrel{~}{R}}_{n}^{\prime 2}}$ is the fraction of the enveloping circle that is occupied by the nine channels. For sufficiently small $\overline{\alpha}$, the nine circular channels and minimal protein spacers (at least 1 nm wide) all fit while touching the DT. In that case: ${R}_{n}^{\prime}={R}_{n}={R}_{dt}+2\overline{\alpha}$. With $R}_{dt$= 8 nm, this is possible up to $\overline{\alpha}\approx$ 3.4 nm. For larger $\overline{\alpha}$, the spacer requirement determines the outer radius of the composite of 9 channels and ${R}_{n}^{\prime}=\overline{\alpha}+\frac{\overline{\alpha}+s/2}{\mathrm{sin}(\pi /9)}$, where $s$=1 nm is the spacer width. ${R}_{dt}$ does not occur in this equation, because the cylindrical channels can no longer (all) touch the DT.
Appendix 6
PDinsight
PDinsight is written in python 3. If available, it uses the numpy module, but does not strictly depend upon that. The program has different modes for computing the parameter requirements for a given $P(\alpha )$ 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(\alpha )$, the mode computeVals will be useful. This computes the expected $P(\alpha )$ values when taking all parameters at face value. Comparison with the subnano 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(\alpha )$. If estimates of PD density and distribution are missing, the mode computeUnitVals, which computes $\mathrm{\Pi}(\alpha )$, could be useful. In this mode, comparison with the subnano channel is impossible.
In tissue level experiments, $P(\alpha )$ is typically measured, but not all ultrastructural parameters will be known. It is likely that PD density $\rho $ and radius (${R}_{n}$, 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 ${R}_{n}$/$\overline{\alpha}$ or $\rho $, 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\ge$ 200 nm), it may be worth exploring the effects of increased central radius (${R}_{c}>{R}_{n}$). This is currently only possible in modes computeVals and computeUnitVals.
Major modes
The core of PDinsight is computing effective permeabilities ($P(\alpha )$) 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(\alpha )$. In mode computeDens, required densities ($\rho $) are computed for given values of maximum particle size $\overline{\alpha}$ (and other parameters). In mode computeAperture, required apertures, given as $\overline{\alpha}$ as well as neck radius ${R}_{n}$, are computed for given values of $\rho $ (and other parameters). In mode computeRnDensityGraph, ${R}_{n},\rho $ curves are computed that together yield a target $P(\alpha )$. The corresponding values of $\overline{\alpha}$ are also reported. These curves can be visualized using any plotting program. Mode sensitivityAnalysis computes socalled elasticities (normalized partial derivatives) around a given set (or sets) of parameters. These elasticities tell how sensitive calculated values of $P(\alpha ),\mathrm{\Pi}(\alpha )$ and constituents like ${f}_{ih}$ are on the parameters involved.
Auxillary modes
In computing $P(\alpha )$, correction factor ${f}_{ih}$ is automatically included. For specific cases such as modelling studies, however, it may be useful to calculate ${f}_{ih}$ separately. For this purpose, several modes exist for exploring inhomogeneity factor ${f}_{ih}$: computeFih_subNano (function of $\overline{\alpha}$; also output values for subnano model), computeFih_pitField_dens (function of $\rho $), computeFih_pitField_xMax (function of $\overline{\alpha}$) and computeTwinning (function of ${\rho}_{pits}$, cluster density).
By default, computations are performed for the unobstructed sleeve model. Most computations can also be performed for the subnano channel model (Liesche and Schulz, 2013; Comtet et al., 2017). Using switch compSubNano, values for the subnano model are also computed.
Graphical user interface
: The GUI to PDinsight is written to facilitate the creation of parameter files and also has a button to run PDinsight directly based on the parameters displayed. In contrast to the generic parameter file, the GUI only shows fields for parameters that are actually used in a specific run mode. The first step of using the gui is to select the run mode (choose from: computeVals, computeUnitVals, computeDens, computeAperture, computeRnDensityGraph and sensitivityAnalysis), followed by ‘load default parameters’. This overwrites any user input from previous modes. All required parameters are shown as text entry fields, with radio buttons for the relevant options. Basic validation occurs on the fly (valid input type, etc). When done, click ‘Run PDinsight’ to write a parameter file and run the program or ‘Export parameter file’ to write a parameter file only. Note that switching modes requires clicking ‘load default parameters’, that is a fresh start.
Additional information can be obtained by clicking the ‘Info’ button and, for certain parameters, clicking on the parameter name. If additional information is available, the mouse cursor changes into a question mark.
Command line usage
(linux): python PDinsight_vXX.py PARAMETERFILE.
(windows): first open a command prompt window (e.g., by searching for ‘cmd’ in the search bar) and go to the directory where PDinsight is located. Then type: py PDinsight_vXX.py PARAMETERFILE
Using the GUI from the command line (linux): python PDinsightGUI_vYY.py or (windows): py PDinsightGUI_vYY.py The correct version of PDinsight_vXX.py should also be available in the current directory.
In the above, XX and YY should be replaced by the respective version numbers.
By default, all outputs are tab separated files (.tsv). The other option is comma separated (.csv), which can be set through the parameter outputType.
The latest version of PDinsight + documentation can be downloaded from https://github.com/eedeinum/PDinsight. DOI: 10.5281/zenodo.3536704
Data availability
PDinsight can be downloaded from GitHub: https://github.com/eedeinum/PDinsight. Documentation on the use of PDinsight.py 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.
References

Emerging models on the regulation of intercellular transport by plasmodesmataassociated calloseJournal of Experimental Botany 69:105–115.https://doi.org/10.1093/jxb/erx337

Tightening the pores to unload the phloemNature Plants 5:561–562.https://doi.org/10.1038/s4147701904508

Plasmodesmata: gateways to local and systemic virus infectionMolecular PlantMicrobe Interactions 23:1403–1412.https://doi.org/10.1094/MPMI05100116

Symplastic intercellular connectivity regulates lateral root patterningDevelopmental Cell 26:136–147.https://doi.org/10.1016/j.devcel.2013.06.010

Homogenization of boundary conditions for surfaces with regular arrays of trapsThe Journal of Chemical Physics 124:036103.https://doi.org/10.1063/1.2161196

On the hydrodynamics of plasmodesmataJournal of Theoretical Biology 74:33–47.https://doi.org/10.1016/00225193(78)902886

Comparison between two methods to measure translational diffusion of a small molecule at subzero temperatureJournal of Agricultural and Food Chemistry 43:2887–2891.https://doi.org/10.1021/jf00059a022

Phloem loading through plasmodesmata: a biophysical analysisPlant Physiology 175:904–915.https://doi.org/10.1104/pp.16.01041

Hindrance factors for diffusion and convection in poresIndustrial & Engineering Chemistry Research 45:6953–6959.https://doi.org/10.1021/ie051387n

Substructure of freezesubstituted plasmodesmataProtoplasma 169:28–41.https://doi.org/10.1007/BF01343367

A general formulation of alternating direction methodsNumerische Mathematik 6:428–453.https://doi.org/10.1007/BF01386093

Intercellular protein movement: deciphering the language of developmentAnnual Review of Cell and Developmental Biology 30:207–233.https://doi.org/10.1146/annurevcellbio100913012915

Photoinducible DRONPAs: a new tool for investigating cellcell connectivityThe Plant Journal 94:751–766.https://doi.org/10.1111/tpj.13918

Plasmodesmata: pathways for protein and ribonucleoprotein signalingThe Plant Cell 14 Suppl:S303–S325.https://doi.org/10.1105/tpc.000778

Modeling the hydrodynamics of phloem sieve platesFrontiers in Plant Science 3:151.https://doi.org/10.3389/fpls.2012.00151

Identification of a developmental transition in plasmodesmatal function during embryogenesis in Arabidopsis thalianaDevelopment 129:1261–1272.

Trapping by clusters of trapsPhysical Review E 61:6302–6307.https://doi.org/10.1103/PhysRevE.61.6302

Diffusion in quasionedimensional structures with a periodic sharp narrowing of the cross sectionRussian Journal of Physical Chemistry B 3:313–319.https://doi.org/10.1134/S1990793109020225

Diffusion in a tube of alternating diameterChemical Physics 370:238–243.https://doi.org/10.1016/j.chemphys.2010.04.012

Plasmodesmata: structure, function and biogenesisCurrent Opinion in Plant Biology 11:680–686.https://doi.org/10.1016/j.pbi.2008.08.002

Plasmodesmata  membrane tunnels with attitudeCurrent Opinion in Plant Biology 14:683–690.https://doi.org/10.1016/j.pbi.2011.07.007

Hindered sedimentation, diffusion, and dispersion coefficients for brownian spheres in circular cylindrical poresJournal of Colloid and Interface Science 124:269–283.https://doi.org/10.1016/00219797(88)903487

Shaping intercellular channels of plasmodesmata: the structuretofunction missing linkJournal of Experimental Botany 69:91–103.https://doi.org/10.1093/jxb/erx225

Symplastic communication in organ formation and tissue patterningCurrent Opinion in Plant Biology 29:21–28.https://doi.org/10.1016/j.pbi.2015.10.007

Regulation of solute flux through plasmodesmata in the root meristemPlant Physiology 155:1817–1826.https://doi.org/10.1104/pp.110.168187

Restricted diffusion through pores with periodic constrictionsAIChE Journal 32:1039–1042.https://doi.org/10.1002/aic.690320615

Plasmodesmata of Brown algaeJournal of Plant Research 128:7–15.https://doi.org/10.1007/s1026501406774

Molecular characterisation of recombinant green fluorescent protein by fluorescence correlation microscopyBiochemical and Biophysical Research Communications 217:21–27.https://doi.org/10.1006/bbrc.1995.2740

Plasmodesmata: gatekeepers for celltocell transport of developmental signals in plantsAnnual Review of Cell and Developmental Biology 16:393–421.https://doi.org/10.1146/annurev.cellbio.16.1.393
Article and author information
Author details
Funding
European Molecular Biology Organization (ASTF 105  2012)
 Eva E Deinum
Nederlandse Organisatie voor Wetenschappelijk Onderzoek
 Bela M Mulder
Engineering and Physical Sciences Research Council (EF/M027740/1)
 Yoselin BenitezAlfonso
Leverhulme Trust (RPG201613)
 Yoselin BenitezAlfonso
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
The authors thank EM Bayer, W Nicolas, M Glavier and L Brocard from the Bordeaux Imaging Centre, Plant Imaging Platform, ((http://www.bic.ubordeaux.fr), for providing the electron micrograph included in Figure 1. Also we thank Philip Kirk (Centre for Plant Sciences, Leeds) for critical review of the python program.
Version history
 Received: June 3, 2019
 Accepted: November 16, 2019
 Accepted Manuscript published: November 22, 2019 (version 1)
 Accepted Manuscript updated: November 25, 2019 (version 2)
 Version of Record published: January 31, 2020 (version 3)
Copyright
© 2019, Deinum et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
Metrics

 2,030
 views

 326
 downloads

 20
 citations
Views, downloads and citations are aggregated across all versions of this paper published by eLife.
Download links
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)
Further reading

 Cell Biology
 Chromosomes and Gene Expression
Heat stress is a major threat to global crop production, and understanding its impact on plant fertility is crucial for developing climateresilient crops. Despite the known negative effects of heat stress on plant reproduction, the underlying molecular mechanisms remain poorly understood. Here, we investigated the impact of elevated temperature on centromere structure and chromosome segregation during meiosis in Arabidopsis thaliana. Consistent with previous studies, heat stress leads to a decline in fertility and micronuclei formation in pollen mother cells. Our results reveal that elevated temperature causes a decrease in the amount of centromeric histone and the kinetochore protein BMF1 at meiotic centromeres with increasing temperature. Furthermore, we show that heat stress increases the duration of meiotic divisions and prolongs the activity of the spindle assembly checkpoint during meiosis I, indicating an impaired efficiency of the kinetochore attachments to spindle microtubules. Our analysis of mutants with reduced levels of centromeric histone suggests that weakened centromeres sensitize plants to elevated temperature, resulting in meiotic defects and reduced fertility even at moderate temperatures. These results indicate that the structure and functionality of meiotic centromeres in Arabidopsis are highly sensitive to heat stress, and suggest that centromeres and kinetochores may represent a critical bottleneck in plant adaptation to increasing temperatures.

 Cell Biology
Highaltitude polycythemia (HAPC) affects individuals living at high altitudes, characterized by increased red blood cells (RBCs) production in response to hypoxic conditions. The exact mechanisms behind HAPC are not fully understood. We utilized a mouse model exposed to hypobaric hypoxia (HH), replicating the environmental conditions experienced at 6000 m above sea level, coupled with in vitro analysis of primary splenic macrophages under 1% O_{2} to investigate these mechanisms. Our findings indicate that HH significantly boosts erythropoiesis, leading to erythrocytosis and splenic changes, including initial contraction to splenomegaly over 14 days. A notable decrease in red pulp macrophages (RPMs) in the spleen, essential for RBCs processing, was observed, correlating with increased iron release and signs of ferroptosis. Prolonged exposure to hypoxia further exacerbated these effects, mirrored in human peripheral blood mononuclear cells. Singlecell sequencing showed a marked reduction in macrophage populations, affecting the spleen’s ability to clear RBCs and contributing to splenomegaly. Our findings suggest splenic ferroptosis contributes to decreased RPMs, affecting erythrophagocytosis and potentially fostering continuous RBCs production in HAPC. These insights could guide the development of targeted therapies for HAPC, emphasizing the importance of splenic macrophages in disease pathology.