From plasmodesma geometry to effective symplasmic permeability through biophysical modelling

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

Obtaining good EM data of PD dimensions is notoriously hard. We therefore opted for a simple geometrical description that allows us to study the effects of PD neck, central region and desmotubule dimensions with as few parameters as possible (see Materials and methods). We modelled a single PD as a 3-part cylindrical channel (Figure 2A), with total length $l$, which would typically equal the local wall thickness. The ends of the channel were modelled by narrow cylinders representing the plasmodesmal ‘neck’ constriction. These have length $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.

Figure 2

Download asset Open asset

We made the arguably simplest choice of modelling particles as (non-additive, 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 so-called steric hindrance reduces the volume that is available for diffusion of the particle’s center in a size dependent manner. Moreover, the maximum particle radius that can pass the PD, $\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: ${\tilde{l}}_n=l_n+\alpha$, ${\tilde{R}}_c=R_c-\alpha$, ${\tilde{R}}_{dt}=R_{dt}+\alpha$ and ${\tilde{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={\tilde{A}}_n/{\tilde{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\widetilde{l_n}\nabla C_n+(l-2\widetilde{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{\tilde{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 so-called 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 slit-pore 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 $\widetilde{A_x}$ with

For completeness, we note that the simplification of a uniform particle flux along the channel axis is violated near the neck-central region transitions, resulting in an error of a few percent (see Materials and methods for further discussion). We now define the permeation constant of a single PD, $\mathrm{\Pi }(\alpha )$, through the rule rule steady-state 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²) 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 3-fold difference when ignoring the particle size dependence of the diffusion constant.

Figure 3 with 1 supplement see all

Download asset Open asset

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.

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 (${\displaystyle 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.

Figure 4 with 2 supplements see all

Download asset Open asset

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{\tilde{l}}_n}$, whereas ${\tau }_{rel}$ ultimately scales quadratically with $R_c$, and, equivalently, linearly with the surface ratio ${\widetilde{\stackrel{~}{A}}}_c/{\widetilde{\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 trade-off 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 ${\tilde{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.

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 ${\displaystyle 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 ${\displaystyle \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 ${\displaystyle 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).

Figure 5 with 1 supplement see all

Download asset Open asset

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.

Figure 6 with 1 supplement see all

Download asset Open asset

Second, we investigated the effect of PDs grouped in small clusters resembling pit fields (see Materials and methods). The average centre-to-centre distance between PDs in pit fields considerably varies across species, with reported/calculated distances between 60 and 180 nm (Terauchi et al., 2015; Schmitz and Kühn, 1982; Danila et al., 2016; Faulkner et al., 2008). The lowest values, however, are from brown algae, which have a different PD structure from higher plants (Terauchi et al., 2012). As a default, we used ${\displaystyle 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 ${\displaystyle p}$ = 5, 6 and 19: indicated with blue lines in Figure 7A; see also Table 1).

Figure 7

Download asset Open asset

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 (${\displaystyle 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).

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

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

PD densities in transverse walls of Arabidopsis thaliana roots were reported by Zhu et al. (1998): vascular: 9.92 ± 0.58, inner cortex: 12.28 ± 0.67, outer cortex: 9.08 ± 0.50 and epidermis 5.42 ± 0.42 PDs/μm². Based on these numbers we assume a PD density of 10–13 PDs/μm² for the tissue level experiment and 5 PDs/μm² 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 ${\displaystyle D}$ = 162 μm²/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 (${\displaystyle 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 ${\displaystyle 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).

Figure 8

Download asset Open asset

Using the model, we also explored the effect of ‘necked’/‘dilated’ PDs by adding a wider central region to PDs. For a central radius ${\displaystyle 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 ${\displaystyle R_c}$ = 26.4 nm) considering a PD density in the order of ${\displaystyle \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₂O₂ treatment. They found a strong decrease in symplasmic permeability to ≈ 1μm/s after treatment with a ‘high’ H₂O₂ 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 (${\displaystyle \rho }$ = 10 μm^–2) or 10.6 nm (${\displaystyle \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₂O₂ 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 ${\displaystyle \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₂O₂ treatment (Figure 8C,D), making DT modification or removal a highly unlikely explanation for this change.

Taken together, these calculations indicate that our model for diffusive symplasmic transport can indeed explain experimentally observed measurements of effective symplasmic permeability, but only with somewhat wider PDs/neck regions than expected yet in line with the observed permeability for GFP and within the range of PD diameters measured in thick cell walls. Alternatively, similar changes in symplasmic permeability can be achieved with several fold higher densities than typically measured. These predictions provide a framework for experimental validation. We also compared the results obtained with our unobstructed sleeve model and the sub-nano channel architecture. Using the sub-nanochannel architecture, much larger PD densities would be required to achieve the same $P(CF)$: roughly twice as large for $\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₂O₂ treatment on effective permeability.

Similar to Smith (1986), we assumed the flux is distributed homogeneously within each cross section along the axis of the channel. This results in a simple mapping to a 1D channel, that is that the average local flux (per unit area of cross section) ∼ 1/available cross section surface. This assumption does not hold close to the transition between neck and central region, that is a sharp transition between narrow and wide cylinders. Numerical simulations showed, however, that the error introduced by the assumption of homogeneous flux turned out to less than 4 percent for ${\displaystyle 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 $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 ${\tilde{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).

For 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 narrow-wide-narrow setup. These calculations, however, critically depend on trapping rates at the narrow-wide transitions. We do not have an expression for these, because the DT takes up the central space of the channel, which, contrary to the case of $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.

To 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 ${\widetilde{\stackrel{~}{A}}}_n$ to accommodate hindrance effects and we replace $A_w$ by $1/\rho$. We also ajust PD length: $\tilde{l}=l+2\alpha$ and $L=L-2\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}{2P_{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 4-fold 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 ${\widetilde{\stackrel{~}{A}}}_n$ by our desmotubule corrected ${\widetilde{\stackrel{~}{A}}}_n$, except in $f(\sigma )$. Additionally, we multiply $f(\sigma )$ by $\xi =({\tilde{R}}_n^2-{\tilde{R}}_{dt}^2)/{\tilde{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.

For 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 sub-nano channel model. In this calculation, a single pit field is modelled as a number of PDs on a triangular (or square) grid with a centre-to-centre distance $d$ between nearest neighbours. We then calculate the pit radius, $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{\tilde{R}}_n^2/{\tilde{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 ${\displaystyle 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 ${\displaystyle R_n}$ = 14 nm and 3–7% with ${\displaystyle 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.

Numbers 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² ($\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 ${\displaystyle \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.

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 ${\widetilde{\stackrel{~}{A}}}_c/{\widetilde{\stackrel{~}{A}}}_n$ of the dilated channel as a rescaling factor for the x-axis. 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$:

Appendix 3—figure 1

Download asset Open asset

For large $R_c$, ${\tau }_{rel}$ becomes proportional to ${\displaystyle {\tilde{R}}_c^2}$. From this we derived a rescaling factor for ${\displaystyle {\tau }_{rel}}$ with ${\tilde{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 x-axis, ${\widetilde{\stackrel{~}{A}}}_c/{\widetilde{\stackrel{~}{A}}}_n$, increases faster for larger particles. The reason is that ${\widetilde{\stackrel{~}{A}}}_n$ decreases relatively faster with particle size than ${\widetilde{\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 ${\displaystyle l_n}$ = 25 nm, this occurs at a wall thickness of ${\displaystyle l\approx }$ 100 nm.

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 (${\displaystyle 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 4—figure 1

Download asset Open asset

Required densities or $\overline{\alpha }$ for the sub-nano 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 ${\displaystyle \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 sub-nano 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\tilde{l}+(1-2\tilde{l})/x)/l$. Similarly, ${\rho }_{gate}$ is computed by assuming a 1 nm thick sub-nano 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 ${\tilde{A}}_n$ is the surface area per cylindrical channel and $\xi =\frac{9{(\overline{\alpha }-\alpha )}^2}{{\tilde{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 ${\displaystyle R_{dt}}$= 8 nm, this is possible up to ${\displaystyle \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}{\sin (\pi /9)}$, where ${\displaystyle 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.

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 sub-nano channel model is possible. In principle, all model parameters must be defined, but missing parameters may be explored using lists of possible values or left at a default (e.g., cell length $L$ and a triangular distribution of PDs), as these have little influence on $P(\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 sub-nano 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 (${\displaystyle 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.

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 so-called 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.

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 sub-nano model), computeFih_pitField_dens (function of $\rho$), computeFih_pitField_xMax (function of $\overline{\alpha }$) and computeTwinning (function of ${\rho }_{pits}$, cluster density).