Robust banded protoxylem pattern formation through microtubule-based directional ROP diffusion restriction

In plant vascular tissue development, different cell wall patterns are formed, offering different mechanical properties optimised for different growth stages. Critical in these patterning processes are Rho of Plants (ROP) proteins, a class of evolutionarily conserved small GTPase proteins responsible for local membrane domain formation in many organisms. While the spotted metaxylem pattern can easily be understood as a result of a Turing-style reaction-diffusion mechanism, it remains an open question how the consistent orientation of evenly spaced bands and spirals as found in protoxylem is achieved. We hypothesise that this orientation results from an interaction between ROPs and an array of transversely oriented cortical microtubules that acts as a directional diffusion barrier. Here, we explore this hypothesis using partial differential equation models with anisotropic ROP diffusion and show that a horizontal microtubule array acting as a vertical diffusion barrier to active ROP can yield a horizontally banded ROP pattern. We then study the underlying mechanism in more detail, finding that it can only orient curved pattern features but not straight lines. This implies that, once formed, banded and spiral patterns cannot be reoriented by this mechanism. Finally, we observe that ROPs and microtubules together only form ultimately static patterns if the interaction is implemented with sufficient biological realism.


Introduction
Plants are able to transport water and nutrients from the ground all the way up to the leaves, potentially more than a hundred meters high, thanks to a highly specialised system of vessels known as the xylem [1]. These xylem vessels are formed 5 by cells that deposit a thick secondary cell wall followed by programmed cell death, leaving a hollow tube [2]. The cell wall reinforcements function to withstand the pressures generated during water transport and may be deposited in intricate patterns depending on the type of xylem. In protoxylem, the 10 secondary wall forms bands or spirals, allowing the vessels to stretch with the surrounding tissue, while metaxylem, formed when longitudinal tissue growth has ceased, tends to be more rigid, with only some well-separated pits for radial transport [3,4,5] (Fig. 1A). 15 The deposition of this secondary cell wall is determined by the position of cortical microtubules on the inside of the membrane that direct cell wall depositing cellulose synthase complexes [7,8,9,10] and secretory vesicles [11]. The microtubule pattern, in turn, depends on the localised activity of 20 ROP (Rho of Plants) proteins that recruit effectors capable of influencing microtubule dynamics, which is best demonstrated for metaxylem development [12,13,14] (Fig. 1A). Since ROP involvement has also been indicated in protoxylem differentiation [15] and the same effectors (MIDD1, Kinesin13A) are ex- 25 pressed during protoxylem development [16,17], patterning of pattern harder to obtain with a reaction-diffusion mechanism. A banded pattern can be obtained on a sufficiently narrow periodic domain, when the circumference of the domain is smaller than the wave length of the pattern as suggested for, e.g., animal tails [31]. However, the protoxylem domain is too large for this 60 effect to occur. Therefore, some other mechanism must impose the orientation of the protoxylem ROP pattern.
The protoxylem pattern has to be properly oriented relative to the cell's growth axis, suggesting that information from the orientation of the cell is transmitted to the ROP pattern. An obvious candidate for the agent of this transmission is the microtubule array itself, since it is well-known for its ability to self-organise into a wide variety of structures, including aligned cortical arrays with a transverse orientation [32,33,34]. This idea is consistent with the observation that, during protoxylem 70 development, microtubules reorient to a transverse array before apparent band formation [35]. In addition, a protein interacting with both microtubules and the plasma membrane can cause microtubules to act as a "molecular fence" that physically re-stricts the movement of active ROP and increased expression of 75 this protein results in more flattened spots in metaxylem [36]. Increasing microtubule stability with taxol treatment or overexpression of MAP70 (a Microtubule-Associated Protein) has a similar effect [13,37]. These observations suggest that, during protoxylem patterning, a transverse microtubule array will 80 pose a barrier to any ROPs moving perpendicular to the array orientation, resulting in anisotropic ROP diffusion (Fig. 1C).
Anisotropic diffusion in reaction-diffusion systems has long been known to influence the shapes of patterns in experimental chemical systems [38] and in models used for texture synthesis 85 [39]. More recently, anisotropic diffusion has been proposed as a general mechanism for orienting stripes formed by Turinglike mechanisms, with horizontal diffusion restriction of a fastdiffusing component (comparable to the inactive form of ROP) resulting in horizontally oriented stripes [40]. Interestingly, in 90 the case of protoxylem patterning we need the opposite: vertical restriction of the active form resulting in horizontally oriented stripes.
Here, we use partial differential equation (PDE) models of ROP-based reaction-diffusion systems with anisotropic diffusion to study active ROP diffusion restriction as a mechanism of ROP pattern orientation in protoxylem development. We also employ a functional decomposition of the diffusion tensor to gain additional mechanistic insight into the orienting power of this mechanism. 100 2. Methods

Modelling approach
We start with two simple reaction-diffusion models for small GTPase-based membrane patterning that can generate coexisting spots, stripes, and gaps [6]. Both models are adaptations of the wave pinning model from [41], with the addition of either protein turnover (WPT model; see Fig. 1D) or negative feedback through GAP activation (WPGAP model) to prevent accumulation of all active ROP into a single cluster (see [6] for a mechanistic explanation). The models assume that active 110 ROP is exclusively membrane bound and inactive ROP is exclusively cytosolic. Consequently, we assume only the diffusion of active ROP will be affected by the microtubule diffusion barriers. Therefore, we will focus our investigation on the diffusion tensor D u of active ROP (of concentration u), assuming the diffusion of inactive ROP (of concentration v), active GAP (G), and inactive GAP (g) to be constant and isotropic. The dimensionless WPT model is given by: where b is the constant activation rate, γ the maximum selfactivation rate, n the hill function exponent of self-activation, ξ the active ROP degradation rate, σ the inactive ROP produc-120 tion rate, D v the inactive ROP diffusion coefficient, and D u the dimensionless diffusion tensor for u. To arrive at the dimensionless forms, all diffusion coefficients and tensors were normalised with the (dimension carrying) diffusion coefficient of unrestricted active ROP (D u,max ).

125
The dimensionless WPGAP model is given by: where active ROP promotes GAP activation at rate c, GAPs are inactivated at constant rate d, and D G and D g are the diffusion coefficients of active and inactive GAP, respectively. In the WP-GAP model, the total amounts of ROP and GAP are conserved, so that the average total ROP concentration T and average to-130 tal GAP concentration T g are constants determined by initial conditions.
We used identical values for the parameters for the reaction parts and the unrestricted diffusion coefficients as in [6]. We have previously performed extensive bifurcation analyses 135 of these models [6] to locate parameter regimes where spontaneous patterning can occur from a homogeneous state (socalled Turing regimes). This analysis directed us to relevant parameter sets that remained useful in the case of anisotropic diffusion. 140

Numerical simulations
We used an Alternating Direction Implicit (ADI) algorithm [42] to solve the reaction-diffusion equations on a two dimensional domain (see Appendix A for details on discretisation schemes). We chose a rectangular domain (190x316 or 60x100 145 dimensionless length units) with periodic boundary conditions in the horizontal direction and zero-flux boundary conditions in the vertical direction, representing the membrane of an elongated cylindrical cell, resembling those found in plant vascular tissue. Smaller square domains (95x95 or 63x63) with fully 150 periodic boundary conditions were used for investigating the orientation mechanism. Unless stated otherwise, simulations were initiated at the homogeneous steady state (see [6]) with a small amount of noise added to each integration pixel. To ensure numerical accuracy, small scale tests were performed using 155 reproducible perturbations as described by [43] to determine the optimal time step and pixel size.

Pattern angle analysis
To characterise the orientation of the simulated patterns, we determined the angle of the patterns ϑ with respect to the horizontal axis at each point on the domain. Since the direction of the pattern essentially runs perpendicular to the concentration gradient, we determined ϑ by taking the angle of the direction perpendicular to the gradient, such that ϑ ∈ − 1 2 π, 1 2 π . This way, ϑ = 0 corresponds to a transverse pattern, while ϑ = − 1 2 π and ϑ = 1 2 π correspond to a longitudinal pattern. An average pattern angleθ over N angles was calculated as follows: where arg is the argument or phase of a complex number and ϑ n is the local angle at each integration pixel.

Homogeneous reduction of active ROP diffusion in vertical direction
In the early stages of protoxylem patterning, before ROP activity starts to affect microtubule density, the density of the microtubule array will be more or less homogeneous. Assuming a transverse array orientation, this homogeneous vertical diffusion barrier can be represented by a spatially homogeneous reduction of active ROP diffusion in the vertical direction, yielding the following dimensionless diffusion tensor: where ψ ∈ (0, 1] is the ratio between diffusion coefficients for restricted and unrestricted diffusion. For isotropic diffusion ψ = 1. Because all off-diagonal elements of D u are zero, the diffusion term simplifies to: Simulation results for various values of ψ are shown in Fig. 2. Increasing the ROP production rate (for the WPT model) or the total amount of ROP (for the WPGAP model) changes the native pattern (i.e., the pattern formed with isotropic diffusion) from spots to stripes to gaps [6], a sequence of patterns often observed in similar models [30]. Vertical diffusion restriction of active ROP flattens spots and gaps and imposes a horizontal 170 orientation on stripes. The stronger the restriction (the smaller ψ), the more the pattern is flattened. Ultimately, this turns the stripes into horizontal bands or, occasionally, spirals and forces spots and gaps to merge, resulting in banded patterns. These results show that vertical diffusion restriction of active ROP can 175 impose a pattern of horizontal bands. Furthermore, to obtain such bands, the native pattern type does not need to be striped.
We note that the observed changes of pattern type require that the diffusion restriction only, or at least predominantly, applies to active ROP. This can be understood from a simple scal-180 ing argument: an equal diffusion reduction for both active and inactive forms is equivalent to independently scaling the x-and y-axis, which will only affect the aspect ratio of the pattern. Similar reasoning suggests that a horizontal pattern can also be obtained by horizontal diffusion restriction of the fast-diffusing 185 component only, which has previously been demonstrated for several other reaction-diffusion systems [40].
Since both WPT and WPGAP models behave very similarly, we will continue our analyses with the WPT model only.

Oblique diffusion restriction promotes spiral formation 190
In reality, the initial microtubule array, and therefore the direction in which active ROP diffusion is restricted, does not have to be completely transverse, but may be somewhat oblique. We define φ ∈ [0, π) as the angle of minimal diffusion restriction with the horizontal axis. This rotates the diffusion tensor resulting in the following diffusion term for active ROP (see Appendix B for derivation): where ψ is again the ratio between the diffusion coefficients of active ROP in the restricted and unrestricted direction. Such oblique diffusion restriction promotes formation of spirals ( boundaries dictate a discrete set of preferred angles that allow connecting spirals without changing the distance between bands (Appendix C). In addition, the zero flux boundaries at the top and bottom of the domain demand that the local gradient there is perpendicular to the boundary. Consequently, spiral 205 angles become slightly steeper than expected based on diffusion angle φ and the periodic boundaries alone.

Orientation mechanism and reorientability
We have seen that the restriction of active ROP diffusion in one direction promotes the formation of bands oriented perpen-210 dicular to that direction. Here, we propose a mechanism for this orientation process and investigate the implications of this mechanism on the reorientability of the pattern. Expansion of existing clusters of high ROP activity can occur when diffusion of active ROP results in ROP activation in neighbouring regions 215 through positive feedback. Conversely, clusters of low ROP activity may expand if diffusion of active ROP from neighbouring regions lowers their activity levels such that they escape positive feedback. This suggests that the speed of expansion is related to the active ROP diffusion coefficient, so that faster dif-220 fusion in a certain direction results in faster expansion in that direction. One may hypothesise that this could lead to a change in the orientation of the pattern. From a geometrical consideration (sketched in Fig. 4A), however, it is to be expected that faster expansion in one direction can only change the orienta-225 tion of curves, but not of straight lines.
A decomposition of the diffusion tensor in components along and orthogonal to the pattern gradient shows that indeed only spots and richly curved patterns can be reoriented by this mechanism, but not straight lines or bands. For this decomposition, we define directions z and w, such that z is oriented along the concentration gradient in descending direction and w is perpendicular to z (Fig. 4C). Then we can rewrite the diffusion term from Eq. 5 in terms of z, w, and the local pattern angle ϑ, defined as the angle between the positive x-axis and the wdirection (see Appendix D): In this case, ϑ ∈ [−π, π). The first term resembles a standard diffusion term in the z-direction. Its diffusion coefficient D zz u depends on ϑ and is therefore constant along straight lines, but larger in the unrestricted than in the restricted direction for cir-230 cular patterns (Fig. 4B). The second term depends on the pattern curvature ∂ϑ/∂w, which is positive for convex clusters, such as spots, and negative for concave clusters, such as gaps and zero for straight lines (Fig. 4B). Because function g(ϑ) is always positive, and ∂u/∂z is always negative, the second term always rep-235 resents a removal of active ROP for spots and an addition for  gaps. Since g(ϑ) is largest in the restricted direction, the second term tends to promote the flattening of circular patterns.
As an illustration, we calculated the magnitude of the complete diffusion terms for some stereotypical patterns generated by the hill function: where we used r = 2x + y for a line and r = x 2 + y 2 for a circle. To turn the circle from a spot into a gap we used 240 u gap = 1 − u spot . The results show that expansion of straight line patterns is dependent solely on diffusion down the gradient, which is the same everywhere along the line (Fig. 4C).
Only curved patterns can expand in other ways, allowing their net orientation to change.

245
These results imply that completely non-curved patterns cannot be reoriented even if the orientation of the diffusion restriction changes, while patterns that contain curves remain open to reorientation. Simulations in which the direction of diffusion restriction is altered after straight bands have been formed 250 confirm this (Fig. 5). Even when the native pattern consists of gaps rather than stripes, a banded pattern was maintained upon returning to isotropic diffusion (Supplemental Figure S.2).
We next investigated which kind of perturbations at the pattern level could reorient a banded pattern to better match the 255 orientation imposed by diffusion restriction. Fully formed horizontal bands were locally stable to small amounts of spatial noise after rotating the diffusion restriction angle (Fig. 6A). Furthermore, they were also stable to moderate local vertical shifts ( Fig. 6B-C), but larger shifts could trigger reorientation in the 260 direction imposed by diffusion restriction (Fig. 6D). This reorientation is not simply an artefact of the severe deformation, since the pattern returned to the horizontal band state after applying the same deformation under vertical diffusion restriction (Fig. 6E). 265 These findings all suggest that diffusion restriction can im- Since ∂u/∂z is always negative and g(ϑ) is always positive, ∂ϑ/∂w determines the sign of the second term from Eq. 7. To emphasise the equal magnitude but opposite effect g(ϑ) · ∂ϑ/∂w has for spots and gaps, opposite signs of this term have been plotted for these two cases. (C) Magnitudes of diffusion components for stereotypical concentration profiles. Shown are: the concentration profile (u) with examples of axes z and w, the diffusion term corresponding to that profile (∇ · D u ∇u), and the two separate terms of the decomposed diffusion term (D zz u (ϑ) · ∂ 2 u/∂z 2 and g(ϑ) · ∂ϑ/∂w · ∂u/∂z).
pose an orientation only on patterns with curves. Since patterns start out fairly rugged from a noisy initial condition, the early stages of pattern formation are expected to be particularly susceptible to reorientation. Indeed, the preferred orientation 270 can already be seen very early during de novo pattern formation (Supplemental Figure S.3). However, fully formed patterns with sufficient curvature can also still be reoriented (Fig. 5D).

ROP-dependent ROP diffusion
Initially, the microtubule array may form an approximately homogeneous barrier to the vertical diffusion of active ROP. Ultimately, however, ROP activity is presumed to result in gaps in the microtubule array, which, in turn, determines the final cell wall pattern (through recruitment of downstream targets). This also means that the diffusion restriction will decrease in these gaps. To study this effect, we made the ratio ψ between unrestricted and restricted active ROP diffusion coefficients a function of active ROP concentration u, such that: To start, we assumed a linear relation between ψ and the dimensionless microtubule density ρ (Fig. 7A): where ψ min is the ratio between the maximally restricted and unrestricted active ROP diffusion components and the maximum of ρ is reached in absence of ROP activity. We take active ROP to reduce microtubule density instantaneously via a hill function: where K ρ is the active ROP concentration at which microtubule density is halved and n ρ is the hill exponent (see Appendix E for parameter choices). Combining Eq. 10 and 11 yields a hill function for the relation between active ROP diffusion and active ROP concentration: Our simulations using these relations for an instantaneous effect 275 of ROP on its own diffusion generally resulted in poorer horizontal band formation (Fig. 7B, Supplemental Figure S.4B). Some simulations got closer to horizontal bands than others, but we could no longer find a broad parameter regime that consistently yielded proper band formation, suggesting that the ro-280 bustness of the orientation mechanism was greatly reduced.

Delayed microtubule density reduction
We saw earlier that developing patterns are particularly susceptible to reorientation, while fully formed banded patterns can be quite resistant to it. These results suggest that a more or less homogeneous diffusion restriction may only be necessary during the development of the ROP pattern, allowing patterning of the microtubule array afterwards, without influencing the ROP pattern. Therefore, we investigated the effect of a delayed response of the microtubule density to ROP activity. Taking a first order approximation for the rate of change towards steady state density ρ * , we get: where α is a rate constant. This means that the diffusion ratio ψ will approach its steady state ψ * according to: Simulations with a homogeneous starting density of ρ = 1 (the maximum, so ψ = ψ min ) and small delays (α = 0.1) yielded patterns comparable to those of the no delay scenario (Fig. 7C,

Positive feedback in microtubule density can prevent travelling waves, stabilising a banded pattern
Thus far, we assumed microtubule density to be merely a 295 reflection of ROP activity, possibly with a delay. However, microtubules have their own complex dynamics that may well have a stabilising influence on the patterning process. Therefore, we will consider a potential positive feedback loop acting on the microtubule density. Microtubule growth occurs only at the tips of existing microtubules. In addition, nucleation of new microtubules occurs predominantly from existing microtubules, with predominantly congruent orientation [44,45,46]. Therefore, we can expect stronger increases in microtubule density at places with higher existing density. The resulting positive feed-305 back loop between microtubule density and microtubule nucleation might serve to stabilise local regions of high microtubule density, thereby preventing travelling wave patterns.
To test this idea, we model this positive feedback with a hill function and let ROP activity promote microtubule degradation, so that dimensionless microtubule density ρ obeys the equation: where γ ρ is the maximum microtubule-dependent growth rate, n ρ is the hill exponent, σ ρ is a constant background increase 310 in microtubule density by nucleation and growth, ξ ρ is a constant linear degradation rate, and ζ is a ROP-dependent degradation rate. Here, ρ is scaled with the density at which positive feedback is half its maximum. In addition, we include a diffusion term with coefficient D ρ to maintain a connection be-315 tween neighbouring integration pixels and prevent artefacts. We take D ρ 1 to keep microtubule "diffusion" well below unrestricted active ROP diffusion. Although microtubules don't really diffuse in their entirety, they tend to spread through growth and nucleation. These processes also favour the existing orien-320 tation [45], which could further promote stable bands. Here, we assume isotropic diffusion as a worst case scenario. Finally, we use α as a tuning parameter for the time scale of the microtubule dynamics relative to the ROP dynamics.
The steady state of the homogeneous system is determined by: which allows for bistability as long as n ρ > 1 (Fig. 8A). In that 325 case, the steady state microtubule density may switch from low to high or the other way around as ROP activity changes.
We use a hill function for ψ to ensure that vertical diffusion is close to its maximum (ψ = 1) for low densities, and close to its minimum (ψ = ψ min ) for high densities (Fig. 8A): where K D is the density at which vertical diffusion is halfway between its minimum and its maximum and n D is a hill exponent.

330
Even when starting from a banded pattern, the positive feedback model without bistability (n ρ = 1) at moderate delay (α = 0.001) yielded travelling waves (Fig. 8B), just like we observed for a direct ROP-dependency with a simple delay. However, when we included bistability (n ρ = 2), the banded pattern re-335 mained static.
However, when simulating de novo patterning of a large domain we still observed non-static patterns at α = 0.001, even with bistability ( Fig. 8C-D). For large delays (α = 0.00001), a banded ROP pattern had sufficient time to develop before sig-340 nificant changes in microtubule density occurred and the banded pattern remained stable. Interestingly, a stable banded pattern seems to form more readily in the gap regime (σ = 0.08) than in the stripe regime (σ = 0.06). We also found some static banded patterns, particularly in the gap regime, for small delays 345 (α = 0.1) that were apparently insufficient to trigger travelling waves. These results show that the precise implementation of microtubule dynamics and the interactions between ROPs and microtubules matters and inclusion of more realistic biological details seems to result in more stable banded patterns.

Directional diffusion restriction can orient ROP patterns into bands and spirals
We have shown that restricting active ROP diffusion in a specific direction is sufficient to change the output from a ROP-355 based Turing-style pattern formation mechanism from patterns without specific orientation into banded or spiral patterns with controlled orientation. Given sufficiently large diffusion restriction, bands and spirals can be generated, regardless of whether the native pattern would consist of spots, stripes, or gaps. The 360 parameters that normally distinguish between these regimes then control only the thickness of the bands. This means that the relatively thick active ROP bands that would be needed to generate the relatively thin microtubule bands observed in protoxylem patterning [10] may actually originate from a Turing mech-365 anism that would generate gaps without diffusion restriction. Additionally, the finding that diffusion restriction can essentially turn spot, stripe, and gap regimes into one large band regime may imply an increased robustness of the protoxylem pattern to changes in ROP activity. By comparison, the spotted 370 metaxylem pattern would more easily be disrupted.
At intermediate levels of diffusion restriction, spots, stripes, and gaps can still be distinguished, but they are flattened horizontally. This finding is in line with the experimental observation that higher expression of a protein responsible for a stronger 375 microtubule barrier effect results in more flattened pits in metaxylem [36]. Therefore, protoxylem patterning may be very similar to metaxylem patterning, with a stronger diffusion barrier effect to obtain a banded pattern and a higher ROP expression to obtain thicker ROP bands. In addition, our finding that the 380 angle of diffusion restriction largely determines the angle of de novo banded pattern formation, predicts that the initial orientation of the microtubule array before apparent band formation determines the final orientation of the banded array. While the precise molecular players and interactions involved in protoxylem patterning remain largely unknown, ROP involvement has been implicated [15]. In addition, we show that different ROP-based Turing mechanisms have a very similar response to a directional reduction in activator diffusion. Therefore, the proposed orientation mechanism seems sufficiently general to apply to different reaction-diffusion systems able to generate coexisting spots, stripes, and gaps, of which many variants exist [6,30,29]. However, it is important that diffusion restriction applies predominantly to the active form, because diffusion restriction of all components in the same way will not 395 yield banded patterns with controlled orientation. The inactive form must therefore be able to bypass the microtubule barriers. This requirement is naturally fulfilled by the predominantly cytosolic localisation of the inactive form.

Fully formed banded patterns are relatively robust struc-
400 tures compared to more curved patterns While anisotropic diffusion has long been known to affect the shape of patterns generated by reaction-diffusion systems [38,39], we not only demonstrated its potential for generating bands and spirals, but also investigated the underlying orienta-405 tion mechanism. Our decomposition analysis of the diffusion term suggests that if ROP cluster expansion is related to diffusion of the active form, any curved structures tend to expand more in the direction in which diffusion is unrestricted. However, straight patterns cannot be reoriented in this way, making 410 fully formed bands and spirals insensitive to any subsequent changes in the direction of diffusion restriction, even in the presence of moderate spatial perturbations.
This robustness of straight patterns has several benefits for the biological system. Firstly, in order to create the actual cell 415 wall pattern, the underlying ROP and microtubule patterns will need to remain approximately stationary for the duration of secondary cell wall synthesis, which may take many hours [10], so the ROP pattern should be robust to any perturbations that may occur in this time frame. Secondly, due to the difficulty of re-420 orienting a pattern of oblique bands, it seems important that the reorientation of the microtubule array to a transverse state occurs before the start of ROP patterning. Experimental observations suggest that this reorientation at least occurs before visible band formation in the microtubule array [35]. Finally, 425 once a transverse microtubule array has contributed to the generation of a banded ROP pattern, the robustness of the banded pattern allows microtubules to disappear in areas of high ROP activity without disrupting the ROP pattern.

430
tion of microtubule dynamics. Our findings show that an instantaneous or very fast response of microtubule density to ROP activity can be disruptive to proper band formation. A sufficiently slow microtubule response, however, gives a robust banded pattern time to form be-435 fore the diffusion restriction is undone. However, this dynamic also results in a kind of delayed negative feedback, which is known for its potential to generate oscillations [47] and travelling waves [48,49] in models. With more realistic microtubule dynamics, we were able to stabilise the banded pattern 440 again. The required difference in time scales of ROP protein and microtubule density dynamics (approximately a factor α) does, however, seem quite large, particularly where the native pattern would consist of stripes. Partly, a large difference may actually be realistic, since changes in microtubule density dur-445 ing protoxylem patterning can take hours [10], whereas GTPase inactivation and removal from the membrane occurs in the order of seconds [50,51]. Also, since the microtubule bands formed in protoxylem patterning are relatively thin [10], the underlying ROP pattern must have relatively thick bands, like those from 450 the gap regime, in which generation of static bands occurred already with smaller delays.
Other biological factors not included here may also help stabilise a banded pattern, potentially reducing the required difference in time scales. Firstly, ROP pattern formation may al-455 ready start before the ROP effectors that interact with microtubules become active, allowing the ROP pattern time to establish without ROP activity having any impact on microtubule density. Transcriptomic data seem to suggest that, in metaxylem, expression of effectors like MIDD1 and Kinesin-13A 460 peaks later than that of ROP11, although the temporal resolution is not very good [52]. Secondly, microtubules are not simply slowly diffusing molecules, but discrete linear structures with complex dynamics and membrane attachment preventing translational movement. Their rigid and linear nature may well 465 contribute to the formation of robust linear bands. Furthermore, microtubule-dependent microtubule nucleation tends to favour the direction of the parent microtubule [45], so that any partial microtubule bands will tend to expand in the direction they already have. Finally, microtubule dynamics are governed largely 470 by processes occurring at their growing or shrinking tips [53]. This suggests ROP activity cannot destabilise existing microtubules except at their ends. Such details of microtubule dynamics have no straightforward implementation in the present modelling paradigm, but detailed stochastic simulation models 475 of the cortical microtubule array with explicit dynamic instability exist [54,55,56,57,58]. Building upon those, a hybrid model combining PDE integration for ROP dynamics with stochastic microtubule dynamics offers an exciting direction for future research on protoxylem patterning. In general, a diffusion term has the following form: where u is the concentration of the diffusing component, D is the diffusion tensor, and D xy = D yx . When the diffusion coefficient is homogeneous (i.e. independent of space), this simplifies to: In this case, the diffusion coefficients do not need to be discretised and discretising the concentrations is straightforward, with u xx at position (x i , y j ) given by: where ∆x is the spatial step size in the x-direction. The discretisation of u yy is done in the same way. For cases with oblique diffusion restriction, u xy can be discretised as follows: For periodic boundary conditions, the last grid point in the periodic direction is simply treated as if it were attached to the first grid point in that direction. For zero-flux boundary conditions, an extra virtual grid point is added at the ends with the same concentration as the grid point before, such that: or the equivalent in the y-direction, where N is the number of the last grid point. As long as at least directions has periodic boundary conditions, this scheme will guarantee mass conservation.

490
For inhomogeneous diffusion, the diffusion coefficients from Eq. A.1 also need to be discretised, because they depend on x and y. For the diffusion term in the x-direction, this was done as follows: where D represents D xx . The diffusion term in the y-direction can be discretised in the same way. This scheme requires the concentrations to be known at the grid points and the diffusion coefficients between grid points (Fig. A.1). To determine the diffusion coefficients between grid points, the average of the 495 two nearest neighbouring grid points was used. For periodic boundary conditions, an extra diffusion node can be added, connecting the last and the first grid points. Zero-flux boundary conditions can be imposed in the same way as before, or by setting the diffusion coefficient between the last grid point and the 500 virtual grid point equal to zero. If oblique diffusion restriction were to be used as well, the xy-terms would also need to be discretised. This can be done as follows: where D represents D xy . The other xy-term can be discretised in the same way. This scheme also requires concentrations at half steps in both directions (Fig. A.2). These can be calculated by taking the average of the four nearest neighbours. For 505 periodic boundary conditions these averages also have to be determined between the last and first grid points. For zero-flux boundary conditions, additional virtual concentration points are required that have the same value as their neighbours. Again, these schemes will guarantee mass conservation as long as at 510 least one of the directions has periodic boundary conditions.

B. Diffusion tensor rotation
For anisotropic diffusion in two dimensions, the diffusion term has the following form: where D xx u , D xy u , D yx u , and D yy u are the components of diffusion tensor D u , with D xy u = D yx u . If diffusion is homogeneous, this simplifies to: The components of the diffusion tensor can be derived from the diffusion rates in the unrestricted and restricted directions (1 and ψ, respectively) and the angle φ of the unrestricted direction with respect to the horizontal direction. For unrestricted diffusion along the x-axis and restricted diffusion along the yaxis, we have: To obtain the diffusion tensor for the case where the unrestricted direction is rotated at an angle φ with respect to the x-axis, we need to rotate D u using rotation tensor R(φ): This way, we can obtain the diffusion tensor D u : (B.5)

C. Preferred spiral angles
With a strictly vertical reduction in the diffusion of active ROP, a fully horizontal banded pattern may form (Fig. 2). The 515 distance between these bands is determined by the dynamics of the system and the zero flux boundaries at the top and bottom, which demand that the pattern starts and ends with either a peak or a valley. This would mean that the preferred number of bands based on the dynamics alone (n bands ) should fall within a quar-520 ter of a band from the observed number of bands. However, the final steady state is also constrained by the initial emergence of the pattern, which depends on the random noise of the initial condition. Therefore, the final pattern may also contain a suboptimal number of bands, which increases the uncertainty on 525 estimates of n bands . The only real ways around this flaw would be to solve the system on an infinitely large domain, which is not feasible, or to average over many simulations with different random seeds for the initial conditions, which is computationally expensive. Therefore, we make a basic estimate of n bands 530 and its uncertainty using the average, minimum and maximum numbers of bands observed in a limited number of repetitions.
Assuming for the moment that we know the preferred number of bands, we can calculate the preferred distance d between the bands as: where H is the height of the domain. Assuming this distance is maintained in a spiral pattern, the angle ϑ of the spiral pattern depends on distance d and the vertical distance between spiral bands d y : Vertical distance d y is constrained by the periodicity of the spiral. The spiral covers a vertical distance D y as it wraps around a domain of width W. This vertical distance D y must be divisible by the vertical band-band distance d y , such that: where n = 1 corresponds to a single spiral, n = 2 to a double spiral, etc. From equations C.1, C.2 and C.3, we can calculate the admissible angles of the spiral pattern if the distance between bands is determined solely by the system dynamics and not by boundary effects: Taking into account the fact that our estimate of n bands may be off by a quarter and the possibility of a suboptimal number of bands, we can determine an expected range for each admissible angle, using the minimum and maximum number of bands observed in simulations with vertical diffusion restriction:

D. Diffusion term decomposition
To study the influence of a pattern's geometry on the effect of vertical diffusion restriction, we decompose the diffusion term in directions z and w. For any point on the domain, direction z runs in opposite direction of the gradient of u, so that derivative u z = − |∇u|. Direction w runs perpendicular to z, so that u w = 0. To decompose the diffusion term, we ultimately need to express the second order spatial derivatives u xx and u yy in terms of z and w. The vectors z and w associated with these directions depend on the pattern angle ϑ, which we define as the angle between positive x-axis and the w-direction: z = sin(ϑ) − cos(ϑ) w = cos(ϑ) sin(ϑ) . (D.1) Using these vectors, we can express derivatives u z and u w in terms of u x and u y : Applying the product and chain rules gives: u xx = sin 2 (ϑ)u zz + sin(ϑ) cos(ϑ)ϑ z u z + cos(ϑ) sin(ϑ)u zw + cos 2 (ϑ)ϑ w u z u yy = cos 2 (ϑ)u zz − cos(ϑ) sin(ϑ)ϑ z u z − cos(ϑ) sin(ϑ)u zw + sin 2 (ϑ)ϑ w u z .

(D.7)
Since z is the direction of the gradient, the pattern angle does not change in this direction, so ϑ z = 0. Also, u zw = u wz = [u w ] z = 0 z = 0. Therefore, the second and third terms of both equations vanish, leaving: u xx = sin 2 (ϑ)u zz + cos 2 (ϑ)ϑ w u z u yy = cos 2 (ϑ)u zz + sin 2 (ϑ)ϑ w u z . (D.8) Substituting these derivatives into Eq. 5 gives: ∇ · D u ∇u = D zz u (ϑ) ∂ϑ ∂w ∂u ∂z D zz u (ϑ) = sin 2 (ϑ) + ψ cos 2 (ϑ) g(ϑ) = cos 2 (ϑ) + ψ sin 2 (ϑ). (D.9) E. Parameters of models with ROP-dependent ROP diffusion 535 Parameters used for the functions linking ROP activity and ROP diffusion are given in Table 1. For the instantaneous ROP effect (Fig. 7B), the maximally restricted diffusion coefficient (ψ min ) was set to 0.25 and other parameters were chosen such that the final value of ψ would be close to its minimum for low 540 ROP activity and close to its maximum for high ROP activity. For the simple delayed variant (Fig. 7C-E), the same value of ψ min was taken and the initial (homogeneous) microtubule density ρ 0 was set such that vertical ROP diffusion would start at its minimum (ψ = ψ min ). For the model with positive feed-545 back in microtubule density (Fig. 8), we tuned the parameters of the density equation to get a bistability in microtuble density that could be switched at low or high ROP levels. We used the same initial microtubule density (ρ 0 = 2) for all simulations. To keep the starting conditions of the vertical active ROP diffu-550 sion coefficient close to that of used for the simple delay, ψ min was lowered to 0.1, with K D = 1. To keep vertical ROP diffusion close to its maximum at low microtubule densities, hill exponent n D was set to 2, creating a sigmoidal curve.