From plasmodesma geometry to effective symplasmic permeability through biophysical modelling

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


Introduction
The formation of spatial patterns in plants requires the transport and interaction of molecular signals. This sharing of information coordinates cell fate decisions over multiple cells and the isolation of cell fate determinants within a cell or group of cells on the same developmental path. Small molecules such as sugars, peptides, hormones and RNAs move long and short distances to coordinate cell/ organ development (Otero et al., 2016). Cell-to-cell transport of proteins, such as transcription factors, is also important in the regulation and/or developmental reprogramming of local cellular domains . A well studied example is SHORT-ROOT (SHR), an Arabidopsis thaliana GRAS family transcription factor, that moves from the stele to cortical-endodermal tissue layers to specify cell fate and root patterning (Nakajima et al., 2001;Spiegelman et al., 2018;Wu and Gallagher, 2013;Wu and Gallagher, 2014). Other mobile factors with developmental importance include TARGET OF MONOPTEROS 7, PEAR transcription factors and miRNAs (Lu et al., 2018;Miyashima et al., 2019;Skopelitis et al., 2018).
Plant cells are connected by channels named plasmodesmata (PDs) that facilitate the transport of these molecules. PD are narrow membrane lined structures embedded in cell walls to allow for symplasmic (cytoplasm-to-cytoplasm) molecular flux (Figure 1). The ER forms a tubular structure called desmotubule (DT) that traverses the PD, leaving a discrete cytosolic sleeve (also called 'cytoplasmic sleeve' in the literature) where molecular transport occurs (Nicolas et al., 2017a;Sager and Lee, 2018). In the region closest to the PD entrances, the cytosolic sleeve appears constricted (neck) in most tissue types, although there are recent observations of 'straight' PDs in meristematic root sections (Nicolas et al., 2017b). Cell walls at PD locations play a key role in regulating its dimensions. The accumulation of callose, a cell wall beta-1,3 glucan polysaccharide synthesized by callose synthases and degraded by bÀ1,3-glucanases (Zavaliev et al., 2011;Amsbury et al., 2017), is the best understood mechanism for the control of PD dimensions and symplasmic transport capacity (i.e. effective symplasmic permeability). Other factors such as membrane composition, shape and number of PDs change during development and between cell types adding extra dimensions to PD regulation (Nicolas et al., 2017a). Mutants blocked in PD form and function are embryo or seedling lethal, highlighting the importance of these structures for normal plant development (Kim et al., 2002;Benitez-Alfonso et al., 2009;Xu et al., 2012).
Small molecules can move via PD by diffusion (non-targeted transport). This is considered to be predominantly symmetrical (Schö nknecht et al., 2008;Maule, 2008), while in certain tissues, such as secreting trichomes (Waigmann and Zambryski, 1995;Gunning and Hughes, 1976) and the phloem (Ross-Elliott et al., 2017;Comtet et al., 2017), hydrodynamic flow may create directionality. The maximum size of molecules that can move by this generic 'passive' pathway is often referred to as the 'size exclusion limit' (SEL), which obviously depends on PD properties and structural features (Dashevskaya et al., 2008). Large molecules can move through PD via an 'active' or 'targeted' pathway overriding the defined SEL. This may involve additional factors that temporarily modify these substrates, target them to the PDs, or induce transient modifications of the PDs to allow for the passage of larger molecules in a highly substrate dependent fashion Maule et al., 2011).
Computational modelling approaches have been applied to model PD transport but, so far, these have mainly focused on hydrodynamic flow and the specific tissues where that matters (Blake, 1978;Bret-Harte and Silk, 1994;Jensen et al., 2012;Ross-Elliott et al., 2017;Comtet et al., 2017;Foster and Miklavcic, 2017;Couvreur et al., 2018). The few existing studies on diffusive transport do not consider neck constrictions or the approach to PDs from the cytoplasmic bulk. Most models consider PDs as straight channels, with advective/diffusive transport through an unobstructed cytosolic sleeve and typically, but not always, account for reduced diffusivity inside these narrow channels compared to the cytosol (Bret-Harte and Silk, 1994;Liesche and Schulz, 2013;Dö lger et al., 2014;Ross-Elliott et al., 2017;Couvreur et al., 2018). Only the oldest of this set, (Blake, 1978), uses a dilated central region in its calculations, but is entirely focused on hydrodynamics. In specific contexts, also a few other geometries are considered. (Ross-Elliott et al., 2017) also consider 'funnel' shaped PDs, which are observed in the phloem unloading zone, but ignore the DT in their diffusion model, as they only calculate a best case scenario for diffusive transport. In the context of size selectivity for small (sugar) molecules in phloem loading, also the so-called 'sub-nano channel model' of PD geometry has been considered (Liesche and Schulz, 2013;Comtet et al., 2017). In this model, symplasmic transport is modelled to be confined to nine cylindrical channels spanning the PD. This was based on a 9-fold rotational symmetry in enhanced 'top view' electron micrographs but never validated experimentally in longitudinal sections. Instead, sparsely spaced axial spoke structures have been reported (Ding et al., 1992;Nicolas et al., 2017b).
Experimental measurement of the parameters that determine effective symplasmic permeability is difficult and many examples exist of misleading and/or conflicting results. Generally speaking two main approaches are used, providing results at different scales that are hard to reconcile. On the one hand, ultrastructural observations using transmission electron microscopy (EM) can provide useful data on PD dimensions and structural features but, despite recent advances, sample preparation affects the integrity and dimensions of PDs to an unknown extent potentially yielding an underestimation of relevant parameters (Nicolas et al., 2017b). On the other hand, tissue level measurement of symplasmic fluxes is achieved using symplasmic molecular reporters, but this is either invasive or limited to few molecular sizes, mostly fluorescein and its chemical relatives (hydrodynamic radii of about 0.4 to 0.6 nm) and GFP derived fluorescent proteins (such as DRONPA-s (28 kDa), Dendra2 (26 kDa), (photoactive and non-photoactive) single GFP (27 kDa, hydrodynamic radius 2.45-2.82 nm) and its multiples [Calvert et al., 2007;Terry et al., 1995;Chudakov et al., 2007;Gerlitz et al., 2018;Kim et al., 2005;Rutschow et al., 2011]). In all cases, the tissue geometry and varying degrees of vacuoloarization can severely complicate the interpretation of the measurements in terms of effective wall permeability for symplasmic transport. Old data on symplasmic permeability use either microinjection or particle bombardment, which allow for a much wider size range of dyes/ molecular reporters, but these techniques can produce cellular stress, which affects PD function (Liesche and Schulz, 2012). Even when using the same dye/fluorescent molecule and the same tissue, these approaches deliver much lower permeabilities than less invasive techniques, demonstrating that they are unreliable for estimating permeabilities in unperturbed plants (e.g. see Haywood et al., 2002, or compare Rutschow et al., 2011and Goodwin et al., 1990. Less invasive methods involve transgenic lines expressing fluorescent proteins under cell-specific promoters (Roberts et al., 2001;Stadler et al., 2005a), which are very time consuming to generate, or photoactivation and photobleaching techniques (Rutschow et al., 2011;Gerlitz et al., 2018). These approaches have yielded valuable insights, but again, both are limited to few proteins/molecular sizes.
In summary, despite recent advances in the development of probes and techniques, effective symplasmic permeability is difficult to assess directly. The fast response of plants to wounding and other stresses, may render part of the ultrastructurally derived parameters less reliable than others, explaining the frequent observation of apparently incompatible results when modelling diffusive symplasmic transport from ultrastructural measurements. In a multi-species analysis correlating photobleaching and electron microscopy results, (Liesche et al., 2019) were unable to find a universal model for matching measurements at the ultrastructural and tissue levels for different interfaces along the phloem loading pathway, illustrating the need for better models. Ideally, we would be able to integrate the results of the experimental approaches at both levels in a model that considers their limitations in order to get more accurate estimates of effective symplasmic permeability and the underlying structural parameters. This brings us to our central question: what do we need to assume about PD size, number, structure, etc. to be able to reproduce tissue level measurements? Moreover, PD geometry changes during development (Roberts et al., 2001;Fitzgibbon et al., 2013), inspiring our second main question: how do distinct features of PD geometry influence transport properties?
Here, we describe the biophysical properties of diffusive symplasmic transport considering detailed PD structural features (such as the DT and the neck region) and the approach from the cytoplasmic molecular bulk towards PDs that are either evenly distributed or clustered into pit fields (Faulkner et al., 2008) (Figure 1). Inside our model PDs, the entire cytosolic sleeve is available for particle diffusion ('unobstructed cytosolic sleeve model'). We investigate how neck/central region, wall thickness, the presence of a DT and PD clustering into pit fields affect transport characteristics for different particle sizes, adding a functional context to some puzzling recent experimental observations. We also apply our framework to compute effective permeabilities for carboxyfluorescein (CF), a fluorescent dye used routinely to measure changes in symplasmic permeability. Comparing calculated and experimentally measured values, we demonstrate that the relatively high effective CF permeabilities observed by Rutschow et al. (2011) can be explained by our model of diffusive 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 'sub-nano channel model' (Liesche and Schulz, 2013) referred to above. Our calculations demonstrate that multi-scale modelling approaches can integrate results from PD structural dimensions and molecular fluxes and reveal conflicts on these determinations. We, therefore, recommend these should be applied systematically when defining effective symplasmic permeability for a particular tissue/molecule and/or biological context. To facilitate this, we share a python program for computing effective permeabilities from PD geometries as a community resource.

Outline of the model
Our aim is to describe the symplasmic transport properties of a cell wall as an effective wall permeability, that is a single number that could be plugged into tissue/organ level models. For this, we split the transport into two parts: the movement through an individual channel representing a PD and the approach to this channel from the cytoplasmic bulk ( Figure 1). This implicitly assumes a homogeneous cytosol. The basic geometrical terminology that we considered in our calculations is introduced in the cartoon PD shown in Figure 1B. An overview of all mathematical symbols is given in Appendix 1.
Obtaining good EM data of PD dimensions is notoriously hard. We therefore opted for a simple geometrical description that allows us to study the effects of PD neck, central region and desmotubule dimensions with as few parameters as possible (see Materials and methods). We modelled a single PD as a 3-part cylindrical channel (Figure 2A), with total length l, which would typically equal the local wall thickness. The ends of the channel were modelled by narrow cylinders representing the plasmodesmal 'neck' constriction. These have length 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 (non-additive, i.e. not interacting among themselves) hard spheres with radius a. 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 a ( 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, a, 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 a to operationalize the SEL concept in a straightforward manner. To avoid confusion, however, we will consistently write a when referring to our model. We introduced rescaled geometrical parameters to account for the reduced available volume in a compact way:l n ¼ l n þ a,R c ¼ R c À a,R dt ¼ R dt þ a andR n ¼ R n À a. 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: a ¼ ðR n À R dt Þ=2. PDs are modelled using multiple cylinders with a total length l, neck (inner) radius R n and neck length l n , central region (inner) radius R c and desmotubule (outer) radius R dt . B,C: Illustration of the impact of steric hindrance and rescaled parameters. The gray areas of the longitudinal (B) and transverse (C) sections cannot be reached by the center of the particle with radius a (steric hindrance). For a concise description of the available volume and cross section area, we use the rescaled lengthsl n ¼ l n þ a,R c ¼ R c À a,R dt ¼ R dt þ a andR n ¼ R n À a. (C) The cross section area available for diffusion on a transverse section was namedÃ, which depends on the particle radius (a).Ã is the area of the white ring in each cross section. The maximum particle size a is illustrated with a dashed circle. For a particle of size a ¼ a,Ã ¼ 0. (D) In practice, particles spend less time diffusing close to the wall than farther away from it (hydrodynamic hindrance). Consequently, the area close to wall contributes less to diffusive transport, as illustrated with purple gradients. These additional hindrance effects are accounted for inÃ. Considering pure diffusion without particle turnover inside the PD, particle flux through the channel is described by qCxyz qt ¼ Dr 2 C xyz , or in steady state: Dr 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, rC x , is inversely proportional to the available surface area A x , so rC c ¼Ã n =Ã c rC n . The total concentration difference over the PD, DC ¼ C l À C 0 is accordingly distributed over the channel: DC ¼ 2l n rC n þ ðl À 2l n ÞrC c . The steady state molar flow rate QðaÞ through each channel is proportional to the entrance cross section: QðaÞ ¼ ÀDÃ n rC n . Solving these equations for rC 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 HðlÞ<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 l. In our case, l ¼ 2a=ðR x À R dt Þ. In the neck region, l ¼ a= a. For the full expression and behaviour of HðlÞ, see Materials and methods.
As HðlÞ already includes the effect of steric hindrance between wall and particle, we can adjust Equation 2 by replacing every instance ofÃ 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, PðaÞ, through the rule rule steady-state flow rate = permeation constant Â concentration difference, yielding

PðaÞ
QðaÞ We also defined t 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ðaÞ, the quantity that can be estimated using tissue level measurements) follows from the definition J ¼ PDC (steadystateflux ¼ permeability Â densityjump): with , the density of PDs per unit wall area (number/ mm 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 PðaÞ (as described in the Materials and methods). As expected, PðaÞ depends on particle size. Two factors underlie this size dependence, which both affect PðaÞ: 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 =a. Figure 3A and ( Figure 3-figure supplement 1) show that hindrance effects have the strongest impact for particle sizes close to the maximum a, whereas the particle diffusion constant always has a large impact Figure 3B. For example, at R n ¼ R c , the 50+ fold difference between a = 0.1 nm and a = 2 nm is reduced to a 3-fold 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 a). We investigated how both the transport volume (using Equation 2) and transport time (t as above) change when the central region is dilated. To compare channels with neck and dilated central region (12 nm ¼ R n 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 t rel (Figure 4). For a more detailed discussion of t rel and its computation, see Materials and methods and Appendix 2.
We then investigated how both Q rel and t rel change with increasing central region radius R c and how this depends on particle radius a and PD length l ( Figure 4). The panels A and B show that 3. Impact of particle size (radius = a) on single pore effective permeability PðaÞ. (A) Dependence of PðaÞ on neck radius (R n ) and a (different line colours, see legend). The diffusion constant D is inversely proportional to particle size (D ¼ d 1 =a). Dashed lines show PðaÞ considering only steric hindrance, solid lines include all hindrance effects. B: Using the same diffusion constant for all particle sizes instead shows that, once particles can pass easily, the particle size dependence of PðaÞ is largely due to the relation between particle size and diffusion constant. Parameters for calculations: l = 200 nm, l n = 25nm, R dt = 8 nm, R c = 17.5 nm. For simplicity we use d 1 = 1 nm 3 /s in this figure. Therefore, only the relative values of the unit permeabilities have meaning (consequently expressed in arbitrary units [a.u.]). The online version of this article includes the following figure supplement(s) for figure 3: 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 (a, 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 l 2ln , whereas t rel ultimately scales quadratically with R c , and, equivalently, linearly with the surface ratioÃ c =Ã n (see Appendix 3 and Appendix 3-figure 1). Figure 4. Impact of central region dilation on molar flow rate (Q) and mean residence time (t). The same legend shown in C applies to all panels. Narrow channels have R n ¼ R c ¼ 12 nm, whereas for necked/dilated channels, R n = 12 nm but R c varies. (A-C) Red curves show the relation between molar flow rate in dilated PD vs narrow PD Q rel ¼ Q dilated ðR n ; R c Þ=Q narrow ðR n Þ whereas cyan curves show the relation between mean residence time in dilated PD vs narrow PD: t rel ¼ t dilated ðR n ; R c Þ=t narrow ðR n Þ. Both quantities are computed for different particle sizes (solid: a » 0, dashed: a = 0.5 nm, sparse dashed: a = 1 nm, dash-dotted: a = 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 t rel has its maximum atl 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 t 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 a. Assuming that control over maximum particle size a 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ðaÞ is proportional to orifice area ( » A n ), we first computed n c ð aÞ, 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 a: Figure 5A displays the n c ð aÞ as a function of the maximum particle size. As an example, when a = 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 a. 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 (a= a). 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 a ( . 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 a ( 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 centre-to-centre distance between PDs in pit fields considerably varies across species, with reported/calculated distances between 60 and 180 nm Figure 5. DT increases the cross section surface area available for transport per channel given a maximum particle radius a. (A) The number of cylindrical channels (n c ) that is required to match the total entrance surface of a single channel with R dt = 8 nm and the same maximum particle radius a. (B) Shows the relative area available for transport (A n ) in relation to relative particle size (a= a) when comparing channels with DT and the equivalent number of cylindrical channels. Total surface area is the same. Solid lines include all hindrance effects (Ã n;dt =Ã n;circle ; cf. Figure 2D). Dashed lines includes steric effects only (Ã n;dt =Ã n;circle ; cf. Figure 2C). The online version of this article includes the following figure supplement(s) for figure 5:  (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 ). 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 Figure 6. Correction factor f ih for inhomogeneous wall permeability depends on PD distribution, cell wall thickness and neck radius. (A) The cartoon shows the geometrical considerations and parameters used to model the diffusion towards PDs. Cell wall inhomogeneity is incorporated as a correction factor f ih , 0<f ih 1, which measures the relative impact of cytoplasmic diffusion towards the locations of the PDs in the cell wall compared to reaching a wall that is weakly but homogeneously permeable (i.e., with f ih ¼ 1). The cytoplasm is considered homogeneous. Each bit of cytoplasm can be assigned to the PD closest to it. With PDs on a regular triangular grid, the cytoplasm belonging to a single PD, with an outer (neck) radius R n , is a hexagonal column with cross section area Aw and 1/2 of the cell length L on either side of the wall. (B-D) f ih is represented as a function of R n . The presence/ absence of DT does not affect the values of f ih ( Figure 6-figure supplement 1A). In all cases, solid lines correspond to: l = 100 nm, L = 10 mm, a = 0.5 nm, a PD density of = 10 PD/mm 2 , and PDs distributed on a triangular grid. Broken lines show the effects of changes in PD density (B), PD length l (C) and PD distribution (D). The online version of this article includes the following figure supplement(s) for figure 6: Figure supplement 1. f ih is not affected by particle size a, presence of DT, or cell length L. 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 . In both cases, we observe the steepest dependency of f ih on for the clusters with the lowest within cluster PD density (pit fields with p = 5, 6 and 19: indicated with blue lines in Figure 7A; see also Table 1).
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ðaÞ 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ðaÞ. 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 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 mm from the quiescent centre). The authors present two experimental setups: a 'tissue level' experiment in which a whole » 50 mm longitudinal section of the root was bleached (estimated effective permeability 6-8.5 mm/s) and a single cell experiment in which a single epidermal cell was bleached (estimated effective permeability 3.3 ± 0.8 mm/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/mm 2 . Based on these numbers we assume a PD density of 10-13 PDs/mm 2 for the tissue level experiment and 5 PDs/mm 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 mm 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 mm -2 for the tissue level experiment and 19 (14 -23) PDs   (Rutschow et al., 2011) with the default model. This table was generated using PDinsight. A: Required density () for a given a and neck radius (R n ). B: Required a and corresponding R n for a given . C: values required to reproduce PðCFÞ = 25 mm/s. Values computed for a 2x, 3x and 4x increase of are also shown. This is done both for a uniform increase of the density (p ¼ 1) and for (repeated) twinning (p>1) from a uniform starting density (indicated in bold). p is the number of PDs per pit.  mm -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 = 10 mm -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 » 1mm/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 ( = 10 mm -2 ) or 10.6 nm ( = 13 mm -2 ), resulting in a = 1.5 nm or 1.3 nm, respectively. The authors also found a strong increase in permeability to » 25 mm/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 = 10 mm -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 would have to increase substantially ( Figure 8B). As an extreme hypothesis, we also calculated the effects of complete DT removal. The increases in PðaÞ 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 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 a = 3.5-4 nm and even larger for smaller a (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 non-targeted, diffusive symplasmic transport based on the dimensions and distribution of PDs and on the size of the mobile particles. For individual PDs, we used a minimal geometrical description that allowed us to extensively investigate the effects of dilation of the central PD region and the implications of a DT at the PD axis on transport properties. Because PDs are narrow, our calculated effective symplasmic permeabilities were heavily affected by molecular hindrance effects. For the effects of PD distribution, we introduced an 'inhomogeneity factor' 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 a 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 Crawford and Zambryski, 2000). This feature, however, can be exploited for the spreading of viruses (Benitez-Alfonso et al., 2010) and other intracellular parasites such as the fungus Magnaporthe oryzae (Kankanala et al., 2007). Functional PDs without DT (and inner diameter of 10-20 nm) have been reported for the brown algae species Dictyota dichotoma (Terauchi et al., 2012). Due to their very high membrane curvature, DT formation requires curvature-inducing proteins (such as reticulons) and a special lipid composition (Tilsner et al., 2011;Grison et al., 2015;Knox et al., 2015). It is likely that performance benefits of the DT offset these costs and disadvantages and it is therefore under evolutionary selection. Additionally, the connection between DT and ER could result in variable degrees of PD occlusion and hence a potential control mechanism for PD accessibility. Park et al. (2019) have started to explore this concept in the context of pressure regulated PD occlusion.
We have also calculated the performance costs (transport rate) and benefits (transport volume per PD) of having distinct central and neck regions. Whereas the transport time scales quadratically with the radius of the central region (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 narrow-wide-narrow structure, and/or (3) the material properties of cellulosic cell walls and PD cell membranes only allow for a distinctly wider central region if the channel is long enough.
In our model, we naturally define the SEL as a, the maximum particle radius that could fit through the model PD, but experimental determination of this value is difficult and often relies on the transport of detectable, typically fluorescent, molecules such as CF. The limited set of suitable molecules, particularly for non-invasive techniques, introduces a large uncertainty in SEL measurements and hence a. Also other biological factors could lead to an underestimation as well as an overestimation of a. For example, in so-called active symplasmic phloem loaders, such as the cucurbits, sucrose moves symplasmically from bundle sheet cells (BSC) to intermediary cells (IC), where it is polymerized into the larger oligomers raffinose and stachyose, that do not diffuse back in detectable amounts (Haritatos et al., 1996;Liesche and Patrick, 2017). Two explanations have been suggested: (1) a discriminating PD SEL at this interface, which prevents the back transport of raffinose and stachyose (Liesche and Schulz, 2013), or (2) open PDs combined with a directional flow which could be sustained by the xylem flow (Comtet et al., 2017). Only the latter could explain the observed amount of sucrose transport (Liesche and Schulz, 2013;Comtet et al., 2017). This example illustrates that the consideration of a symplasmic flow could largely affect calculated permeabilities and fluxes.
An overestimation of a could occur for non-spherical molecules or temporal variations in PD properties. Although a molecule's hydrodynamic radius is a better predictor of its symplasmic transport efficiency than its molecular weight (Terry and Robards, 1987;Dashevskaya et al., 2008), it conceptually assumes a static replacement sphere. Molecules may be more flexible and/or have a shortest dimension than what is captured by its diffusive behaviour in bulk. PDs might also accommodate molecules that are larger than expected, either through interactions with specific PD proteins (Benitez-Alfonso et al., 2010)  The framework we have developed for so-called 'simple' PDs also provides an intuition for the functional implications of complex geometries such as 'twinned ', 'branched' or 'funnel' PDs (Ehlers and Kollmann, 2001;Ehlers and van Bel, 2010;Faulkner et al., 2008;Ross-Elliott et al., 2017). All else remaining equal, 'twinned' PDs have twice the entrance surface area, which would result in doubling the effective permeability PðaÞ. This increase, however, will be reduced because of the less uniform PD spacing in a density dependent manner ( Figure 7A). 'Branched' or 'complex' PDs contain multiple sub-channels (branches) on at least one side with typically a single shared central cavity connecting all branches (Oparka et al., 1999;Roberts et al., 2001;Fitzgibbon et al., 2013). In the leaf sink/source transition, massive branching is observed and, coincidentally, the number of PDs is reduced (Roberts et al., 2001). The formation of many channels per PD could help to maintain sufficient transport capacity for smaller molecules. If so, the increase in the number of typically narrower channels should be much larger than the decrease in total (simple or complex) PD number. Our computations of 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 (Ross-Elliott et al., 2017) and show a wide opening on the PSE (protophloem sieve element) side and a narrow opening on the PPP (phloem pole pericycle) side. (Ross-Elliott et al., 2017) model these as a triangular funnel that reaches its narrowest diameter only at the (PPP) bottom. There appears to be, however, a longer neck-like region at the narrow end of variable length. As hindrance is by far the highest in the narrowest section, the length of this narrow part would be a vital parameter in correctly estimating the transport permeabilities of these PDs.
We have applied our model to calculate the effective permeability for fluorescein in transverse walls of Arabidopsis root tip cells (Rutschow et al., 2011). Assuming purely diffusive transport and parameters based on various ultrastructural measurements, we were able to reproduce the observed effective permeabilities for CF and to assess the plausibility of different hypotheses aimed at resolving the conundrum of apparently incompatible measurements at different scales. For resolving this conundrum, we assumed that not all PD dimensions are reliably measured with EM. We could reproduce the measured values with somewhat wider PDs/neck regions or several fold higher PD densities than usually measured by EM. Of these, the increased radius seems the more plausible scenario, in line with the requirements for efficient GFP transport reported to occur among root meristem cells (Benitez-Alfonso et al., 2009;Benitez-Alfonso et al., 2013;Nicolas et al., 2017b), and similar to 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 » 2-fold difference in PD density (Zhu et al., 1998). The obvious next step would be testing more data sets of different interfaces/plant species where purely diffusive symplasmic transport is expected. First of all, it would be ideal to test if a near or complete match between tissue level and ultrastructural measurements can be produced if all measurements are performed on the same system with the same growth/treatment conditions. Additionally, more testing could yield a better understanding of potential systematic side-effects of modern EM preparation techniques and/or uncertainties in the tissue level measurements, which would show as systematic vs random required adjustments of the model parameters. A very exciting outcome would be the discovery of distinct clusters in required parameter adjustments that could be related to cell wall properties, PD or interface type, etc. Additional model testing would become easier if the results of tissue level experiments are reported in the form of effective symplasmic wall permeabilities (in mm/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/mm 2 , a reduction from a = 4.2-5.2 nm to a = 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 a = 8.8-10.5 nm, or a ±4-fold increase in PD density (possibly through 2 rounds of twinning/duplication), or less extreme changes if both parameters increase simultaneously (see Table 2C). The required increase in PD density should occur relatively fast, that is within the applied incubation period of 2 hr, and is so large that it should be readily detectable with EM.
The fact is that to reproduce experimentally measured CF effective permeabilities with our model, we had to deviate from ultrastructural based values for at least one parameter. Potential sources for these variations are: (1) ultrastructural studies might underestimate R n because plants could respond to pre-EM manipulation by closing PDs, like they do in response to microinjection or particle bombardment (Haywood et al., 2002;Liesche and Schulz, 2012), (2) PD integrity could be affected during processing for TEM leading to an underestimation of PD densities, (3) the mechanical properties of cell walls and membranes provide a flexibility in the channel that could to some degree accommodate molecules larger than the apparent R n (Abou-Saleh et al., 2018;Yan et al., 2019;Amsbury and Benitez-Alfonso, 2019). For a passive transport mechanism, the elastic energy required for these reversible deformations would have to be in the order of a few 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 (Bret-Harte and Silk, 1994;Liesche and Schulz, 2013;Dö lger et al., 2014;Ross-Elliott et al., 2017). Using typical PD dimensions and no clustering, inhomogeneity factor 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ðaÞ (Appendix 2), but results in parameters that are directly measurable on EM images.
Comparing the unobstructed sleeve architecture to the sub-nano channel architecture, we found that the latter requires roughly twice as high PD densities to produce the same permeability values PðCFÞ in the (Rutschow et al., 2011) experiments. This difference is due to the increased hindrance effects in cylindrical channels vs annular channels with the same cross sectional area. In the future, sleeve models could be refined with the consideration of central spokes (Ding et al., 1992;Nicolas et al., 2017b) and variability of PD dimensions within a single cell wall (Nicolas et al., 2017b;Yan et al., 2019). Simple considerations of the available volume suggest that the addition of spokes will increase hindrance effects, but most likely to a lesser extend than the sub-nano channel structure. Detailed molecular simulations could be a valuable tool to assess this effect.
Other future applications could be the coupling of our detailed PD level calculations of effective symplasmic permeability with tissue level models, which would allow for investigating the impact of microscopic changes on developmental and physiological processes (for example see Foster and Miklavcic, 2017;Couvreur et al., 2018). Depending on the context, it would then be useful or even required to also implement hydrodynamic flow through the PDs. Many ingredients are available for doing this while maintaining the distinguishing features of our mode, including hindrance factors (Dechadilok and Deen, 2006), but as far as we know, the theoretical and numerical results that we use for calculating 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 Similar to Smith (1986), we assumed the flux is distributed homogeneously within each cross section along the axis of the channel. This results in a simple mapping to a 1D channel, that is that the average local flux (per unit area of cross section)~1/available cross section surface. This assumption does not hold close to the transition between neck and central region, that is a sharp transition between narrow and wide cylinders. Numerical simulations showed, however, that the error introduced by the assumption of homogeneous flux turned out to less than 4 percent for l = 200 nm, the shortest l with experimentally observed neck region in Nicolas et al. (2017b) (Appendix 2-figure 1) and will be less for longer channels. This error can be considered irrelevant given the quality of available data on PD dimensions and the many molecular aspects of PD functioning that are necessarily neglected in a simple model.

Hindrance factors
Hindrance factors HðlÞ 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.

Relative molar flow rate and MRT
For assessing the impact of the neck constriction on PD transport, we defined two relative quantities: Q rel ¼ Q dilated =Q narrow and t rel ¼ t dilated =t narrow (Figure 4, Appendix 3-figure 1). Using Equation 2 for QðaÞ, Q rel is well defined: For t rel we first needed an expression for t itself. Ideally, this would be a MFPT, which could calculated in a way similar to t k 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ðaÞ and defined t 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 t rel , the concentration difference cancels from the equation, solving the problem: This method of computing t 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
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) t k through a simplified PD and a column of cytoplasm surrounding it. We then convert t k 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= 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 s, the fraction of the wall occupied by the circular PDs, increases from 0 to 1. f ðsÞ 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 k n must satisfy A w k w ¼ A n k 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 a. 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Ã n to accommodate hindrance effects and we replace A w by 1=. We also ajust PD length:l ¼ l þ 2a and L ¼ L À 2a. At the same time, we adjust f ðsÞ 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 k w : where H c ðlÞ is the hindrance factor for cylindrical pores (see Materials and methods). In the same fashion, we also adjust k n .
We then invert the relation t k ¼ L 2 2D þ L 2Peff , where we write P eff for the effective wall permeability (Makhnovskii et al., 2009), to obtain P eff ¼ L 2t k ÀL 2 =D . With this, we can compute f ih ¼ P eff =ðPðaÞÞ, where PðaÞ 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 t k 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Ã n by our desmotubule correctedÃ n , except in f ðsÞ. Additionally, we multiply f ðsÞ by ¼ ðR 2 n ÀR 2 dt Þ=R 2 n . Numerical calculations in a simple trapping setup confirm the validity of reducing f ðsÞ proportional to the area occupied by the desmotubule whilst calculating s 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 k 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 ðsÞ for different regular trap distributions as given in Berezhkovskii et al. (2006). As A w is calculated implicitly from 1=, no other adjustments were necessary.
Correction factor f ih for pit fields 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 centreto-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 t k 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 k w;pit and t k;pit : where p is the number of PDs per pit and ¼ pR 2 n =R 2 pit is the fraction of available pit area that is occupied by available PD area, and In these equations, is the total PD density. In our graphs, we either keep constant while increasing p to investigate the effect of clustering, resulting in a pit density pits of =p, or keep pits constant to investigate the effect of (repeated) PD twinning. As a default, we used d = 120 nm based on distances measured from pictures in Faulkner et al. (2008) of basal cell walls of Nicotinia tabacum leaf trichomes. To verify our calculations, we compared them with a single step calculation with large circles only, that is with radius R pit and density =p. As results in 3D suggest that for strongly absorbing clusters, the outer radius and cluster density dominate the diffusion (survival time) process (Makhnovskii et al., 2000), this should produce a lower bound to 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 a with default model Numbers in Table 2 are computed based on forward computation of PðaÞ given , a, corresponding R n and other parameters with increments of 0.1 PD/mm 2 () or 0.01 nm ( a etc.) and linear interpolation between the two values that closest match the target PðaÞ. This yields an error of less than 0.0001 mm/s on PðaÞ. We use a = 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. .txt), figure 7 (fig7_pars.txt, fig7d_pars.txt), figure 8b (fig8b_pars.txt) and cross validation of figure 6 (fig6_pars.txt).

Additional information
. Transparent reporting form

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.  Table 2) and the sub-nano channel model. A: Required density () given a and corresponding neck radius (R n ). B: Required a and corresponding R n given . For P (CF) = 25m/s, also values for a 2x, 3x and 4x increase of are computed. This is done both for a uniform increase of the density (p ¼ 1) and for (repeated) twinning (p>1) from a uniform starting density (indicated in bold). A, B: The + sign at R n indicates that the stated R n is too narrow to fit nine sub-nano channels that touch a DT with R dt ¼ 8 nm. Models used for calculating required densities: , R n : default unobstruced sleeve model, sn , R n;sn : sub-nano channel model, sn;neck : sub-nano channel model restricted to the neck regions (l n = 25 nm), and gate : 1 nm thick structures at both PD entrances locally similar to the sub-nano channel model. *: a sn is calculated using R n that allows for 1 nm spacers. C: p is the number of PDs per pit. This table was generated using PDinsight.

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ðaÞ 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ðaÞ, the mode computeVals will be useful. This computes the expected PðaÞ 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ðaÞ. If estimates of PD density and distribution are missing, the mode computeUnitVals, which computes PðaÞ, could be useful. In this mode, comparison with the sub-nano channel is impossible.
In tissue level experiments, PðaÞ is typically measured, but not all ultrastructural parameters will be known. It is likely that PD density 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 / a or , respectively. Additional poorly known parameters can be explored as suggested above. If uncertain, it is strongly advised to explore PD length l ( » cell wall thickness). For thick cell walls (l ! 200 nm), it may be worth exploring the effects of increased central radius (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ðaÞ) 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ðaÞ. In mode computeDens, required densities () are computed for given values of maximum particle size a (and other parameters). In mode computeAperture, required apertures, given as a as well as neck radius R n , are computed for given values of (and other parameters). In mode computeRnDensityGraph, R n ; curves are computed that together yield a target PðaÞ. The corresponding values of a 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ðaÞ; PðaÞ and constituents like f ih are on the parameters involved.

Auxillary modes
In computing PðaÞ, 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 a; also output values for sub-nano model), computeFih_pitField_dens (function of ), computeFih_pitField_xMax (function of a) and computeTwinning (function of pits , cluster density).
By default, computations are performed for the unobstructed sleeve model. Most computations can also be performed for the sub-nano channel model (Liesche and Schulz, 2013;Comtet et al., 2017). Using switch compSubNano, values for the sub-nano model are also computed.