Turing-like mechanism in a stochastic reaction-diffusion model recreates three dimensional vascular patterning of plant stems

Vascular tissue in plants provides a resource distribution network for water and nutrients that exhibits remarkable diversity in patterning among different species. In many succulent plants, the vascular network includes longitudinally-oriented supplemental vascular bundles (SVBs) in the central core of the succulent stems and roots in addition to the more typical zone of vascular tissue development (vascular cambium) in a cylinder at the periphery of the succulent organ. Plant SVBs evolved in over 38 plant families often in tandem with evolutionary increases in stem and root parenchyma storage tissue, so it is of interest to understand the evolutionary-developmental processes responsible for their recurrent evolution and patterning. Previous mathematical models have successfully recreated the two-dimensional vascular patterns in stem and root cross sections, but such models have yet to recreate three-dimensional vascular patterning. Here, a stochastic reaction-diffusion model of plant vascular bundle patterning is developed in an effort to highlight a potential mechanism of three dimensional patterning–Turing pattern formation coupled with longitudinal efflux of a regulatory molecule. A relatively simple model of four or five molecules recreated empirical SVB patterns and many other common vascular arrangements. SVBs failed to develop below a threshold width of parenchymatous tissues, suggesting a mechanism of evolutionary character loss due to changes in the spatial context in which development takes place. Altered diffusion rates of the modeled activator and substrate molecules affected the number and size of the simulated SVBs. This work provides a first mathematical model employing a stochastic Turing-type mechanism that recreates three dimensional vascular patterns seen in plant stems. The model offers predictions that can be tested using molecular-genetic approaches. Evolutionary-developmental ramifications concerning evolution of diffusion rates, organ size and geometry are discussed.


Introduction
The diversity of biological patterning is astonishing. Underlying much of that diversity are conserved genetic modules whose functioning, number, and timing, rate, and location of activity vary from lineage to lineage [1]. Vascular tissues in plants present one suite of biological patterns that function as resource distribution networks that carry water and nutrients throughout the plant body and simultaneously provide mechanical support [2]. The evolution of vascular tissues represents one of the most important adaptations of land plants. Vascular tissues' functions of conduction and support permitted the invasion of otherwise inaccessible arid, terrestrial regions and contributed to vascular plant long distance dispersal and competitive dominance [2,3].
Interest in modeling vascular development has grown recently, and multiple models stemming from a reaction-diffusion framework or molecular pathway analysis have recreated the various arrangements of vascular bundles in plant stems and roots. The classic pattern consists of a cylinder of vascular tissue or radially-arranged vascular bundles near the periphery of stems and roots [4][5][6][7][8][9][10][11][12] with water conducting tissue (xylem) developing proximally to the zone of vascular differentiation (the vascular cambium) and carbohydrate-transporting tissue (phloem) developing distally.
Thus far, mathematical models have not recreated three dimensional patterning that is essential for resource distribution, and many of these models do not explore patterning due to variant cambial activities that alter vascular patterning from the classic pattern. Here, a model is developed to illustrate potential mechanisms of evolution and development of both classic vascular arrangements in stems as well as vascular patterning of supplemental vascular bundles (SVBs) in ground tissue of stems and roots.
Previous modeling work presents common themes that underlie patterning processes in biological systems. In animals, models provided insights about mechanisms of developmental patterning of hair [13], reptile scales [14], feathers [15], and spotting and striping in multiple taxa [16][17][18][19][20]. In plants, mechanisms responsible for the arrangement of leaves on stems (phyllotaxis) [21,22] and the patterning of hairs (trichomes) on leaves [23,24] are also understood at some level. Although it is unlikely that a single mechanism is responsible for the patterning of such repeating, serially homologous units, these models consist of common components: reactions among the molecules, movement of molecules, and boundaries that define limits to the movement of the molecules. Often such models include an activator molecule that initiates development of a serially homologous feature (e.g., hair, trichome, leaf primordium) when it reaches a sufficiently high concentration. Many of the models also include an inhibitor molecule whose production is fostered by the activator and which suppresses both the activator and the development of the feature in adjacent regions.
In 1952, Alan Turing [25] encapsulated these core dynamics in the reaction-diffusion framework. This framework models both the movement of molecules as well as their reaction dynamics using (mathematically) continuous partial differential equations (PDEs). Munteanu and Solé [26] refer to this continuous approach of reaction-diffusion modelling as "the standard tool for evo-devo". Reaction-diffusion models cast evo-devo models in a predictive framework [15], as the developmental, and hence phenotypic, impacts due to evolutionary changes in gene expression rates, molecular degradation rates, diffusion rates, or reaction kinetics can be predicted directly from an analysis (or simulation) of the reaction-diffusion models.
Patterning of supplemental vascular bundles (SVBs) raises interesting developmental and evolutionary questions. SVBs evolved in at least 38 plant families (Fig 1) [27][28][29][30][31] and occur in both stem and root tubers [30]. Like other vascular bundles, SVBs have vascular cambia that produce xylem and phloem in various configurations, but SVBs are separate from the secondary vascular cambium that occurs in a cylinder at the periphery of stems and roots of woodforming plants. In contrast, SVBs tend to occur in the centers of succulent stems and roots or in wide stem cortices of plants such as cacti (Fig 1). They tend to evolve in tandem with evolutionary increases in parenchyma storage tissues [30] and are implicated in distributing resources within these putatively diffusion-limited storage structures [28,30,33]. As such, regions of proliferating storage parenchyma and SVBs are often associated with xeric succulent life forms such as cacti, as well as with plants whose tubers are of agronomic importance such as turnips, kohlrabi, and sweet potatoes.
Within parenchymatous roots and shoots, SVBs are dispersed within the cortices, pith, or, more generally, centers of stems and roots, forming a spotted pattern in cross section and a matrix of (usually) longitudinally-oriented vessels in longitudinal section. In particular, four features of this network of SVBs are distinctive. First, these SVBs are often localized in the pith with zones of parenchyma separating them from the secondary vascular cambium (Fig 2 region B). Second, they tend to have a fairly uniform pattern of dispersion but deviate from complete uniformity. Third, they appear to be relatively rare in thin stems and roots, but more frequent in wide succulent stems, hypocotyls, or roots such as those in Adenia [30,34,35], Cactaceae [32,33], and Brassica (SVBs in turnip roots and kohlrabi stems illustrated in Figs 1 and 2 are absent from thin-stemmed kale and bak choi [36]), and they show a pattern of repeated evolutionary gain and loss in some clades [30]. Fourth, the number, size, density, and spatial distribution of regions differentiating into SVBs vary from taxon to taxon (compare turnip to sweet potato in Fig 1).
A stochastic, reaction-diffusion model consisting of interacting molecules (labeled H, B, M, P, and A) was developed to highlight potential mechanisms of three dimensional vascular patterning. To address the hypothesis that changes in stem and root size alone can cause changes in developmental patterning of SVBs, simulations investigated varying widths of stems. SVBs were expected to develop in wide stems and expected not to develop in narrow stems even when the molecular networks driving pattern formation are identical in wide and narrow stems. To investigate the role that changes in molecular diffusion rates have on patterning, diffusion rates of a regulatory molecule (labeled H) and its associated substrate (B) were varied while all other model parameters were kept constant. Application of Wei and Winter's theory [37] predicted the number and size of SVBs in 2D sections. Their work predicted that the ratio of the diffusion rates of substrate and regulatory molecules most directly influence the density of SVBs, whereas the magnitude of the diffusion rate of the regulatory molecule influences their spatial scale. By altering additional parameters in the model that is developed here, other vascular patterns from a diversity of vascular plant lineages were recreated. Results of these simulations are discussed in terms of potential evolutionary-developmental mechanisms and their connections to known molecular regulators of vascular development.

Model creation
The core 2D model. To derive a reaction-diffusion model of vascular bundle patterning, a "core" model of two reacting components formed a foundation that was augmented, molecule by molecule, to arrive at a four-molecule HBPM model. Ultimately, this model needed additional augmentation with a fifth molecule, A, to correspond better with previously observed biological processes. In this methods section, outcomes of the simpler models and their limitations are described in order to motivate the more complex models.
Gray and Scott [38] and Schnakenberg [39] introduced a set of reaction diffusion models that form the core of the current model. These previous models are recognized for the diversity of spotted and striped patterns that they produce in 2D and 3D spaces. Both the Gray-Scott and Schnakenberg models have the same core activator-substrate dynamics with two interacting components. In the model, a dimer, here labeled H, autocatalyzes its formation in the presence of a 'substrate' (labeled B) that is consumed when it interacts with H. (The [H] 2 term in Eqs 1 and 2 below is interpretable as a bi-molecular, or dimeric, interaction.). A Gray-Scott-Schnakenberg system of equations describing activator-substrate temporal and spatial dynamics is: The k i values are rate constants, the D i values are diffusion rates, and the r 2 Laplacian diffusion operator determines rates of movement between (spatial) regions. Eqs 1 and 2 model the concentrations of H and B as a deterministic system of non-negative, real-valued numbers that can take arbitrarily small positive values. When considering molecules interacting in cells, this framework may not be appropriate, as the number of molecules is whole number-valued, and stochasticity is an intrinsic part of molecular dynamics in biological systems [40], especially when low numbers of molecules (seen in many regulatory molecules) are present in the cell. For example, under some parameterizations, the Schnakenberg model creates 'stochastic resonance' in which spurts of chemical concentrations occur when modeled as a stochastic process, whereas the deterministic model admits a constant steady state; additionally, reaction trajectories can shift between locally stable states in stochastic systems (e.g., the Schlogl model [41]). These stochastic outcomes are implemented as examples in the software developed herein to model vascular evo-devo (S1 Program).
To accommodate the possibility of stochastic effects, the "core" HB model (Eqs 1 and 2) developed here tracked each molecule as it diffused and reacted stochastically according to the above chemical reactions, a departure from previous models of vascular development that used continuous partial differential equations. To implement the stochastic model, the PDE of Eqs 1 and 2 was represented as a set of chemical reactions (Reactions X 1 -X 5 of Table 1). This framework using chemical reaction formulae to model chemical dynamics was illustrated in Turing's 1952 paper [25] and further developed by Gillespie [42]. An exact stochastic simulation using the stochastic simulation algorithm (SSA) of Gillespie [42] produces a realization of the stochastic master equation that models the chemical reactions between B and H.
Relatively recent advances in Monte Carlo sampling make this simulation approach feasible for a large number of reactions [43,44]. S1 Appendix provides details concerning algorithms and implementation. As S1 Appendix describes, all simulations took place in a discretized lattice of equal-volume cubical regions, hereafter called "voxels". This finite element approach is a standard approach to model spatial processes, and a similar finite element approach was used previously to model plant vascular patterning (e.g., [6]). Simulations were carried out in either a circle in 2D (representing a stem cross section) or in a right circular cylindrical solid in 3D (representing a plant stem), and boundaries of the circle or cylinder were reflective (zero-flux Neumann boundary conditions). Reactants were initialized to concentrations of 0.
The reaction-diffusion SSA was implemented in a Java 7 program written specifically for these simulations. This program can run arbitrarily configured stochastic reaction-diffusion systems (S1 Program ZIP archive file includes the application program, source code, examples, and documentation). The program provides a graphical user interface developed using Net-Beans 7.4 IDE and Java Swing (Oracle) to define the reaction-diffusion system, initialize it, and run the SSA. Capabilities are in place to save data and images of the system states, and Quick-Time-formatted movies of the simulation can be produced. The model simulation state can be saved, so the simulation can be stopped and continued at any point in time. As a Java program, it can be run on any platform that supports Java 7 or later.
Model with 3D structure. The two-molecule HB model produced spots of high [H] that were regularly arrayed within the simulated 2D stem cross section, as expected (see, e.g., [37]). In these simulations, these regions of high [H] represent regions of vascular tissue differentiation. Running the deterministic Gray-Scott-Schnakenberg reactions in 3D in the software package Ready v. 0.8 [45] resulted in a connected network of conduits with high [H] (Fig 3A; S1 Table, Simulation A; S1 Model), whereas stochastic implementations of Eqs 1 and 2 resulted in disconnected regions of high [H] in 3D (S1 Movie provides one realization of the 3D implementation of the HB model). Neither the stochastic nor the deterministic approaches reproduced the longitudinal orientation of vascular bundles seen in plant stems.
To improve the model of longitudinal development of vascular bundles, the HB model was augmented to include an additional molecule-an efflux carrier molecule, labeled P. The HBP model is represented by Reactions X 1 -X 8 of Table 1 in which P transported H longitudinally (reaction X 8 : P"#H). H, in turn, promoted expression of P (reaction X 6 : H!H+P). All parameterizations of the 3D simulations presented here that included P are in S1 Table and S2 Table. Molecular transport via carrier molecule P can be shown to follow Michaelis-Menten kinetics as follows. During passive diffusion, a molecule that diffuses outside of its voxel is equally likely to diffuse into any of the 6 adjacent voxels. Two of these six are in the longitudinal orientation, so a molecule that is passively diffusing out of its voxel will diffuse longitudinally with a probability of 1/3. When P is included in the model, the probability that a molecule of H that is leaving its voxel, i, moves longitudinally ("#H) was modeled as (Eq 3): Here, α can be thought of as H's affinity for P and D H is the coefficient of passive diffusion for H. Setting the Michaelis-Menten parameter V max to 2/3 and the Michaelis-Menten parameter K m to 6D H /α, the probability that H moves longitudinally is seen to obey Michaelis-Menten kinetics above the background passive longitudinal diffusion probability of 1/3 by comparison to the Michaelis-Menten equation.
Refining spatial pattern through inhibition of H by M. The augmented model in which P polarized H longitudinally produced longitudinally oriented bundles (Fig 3B-3D; S2 Movie) throughout the stem. In eudicot stems, however, medullary SVBs tend to occur in the centers of stems with a zone of separation between them and the secondary vascular cambium (Fig 2  region B). The three-molecule model (HBP) was further augmented to better match observed distributions of SVBs by including an inhibitor molecule, M. This model that includes four molecules H, B, P, and M (hereafter called the HBPM model; Fig 4) is represented by Reactions X 1 -X 12 of Table 1. M was produced in the outer 4% of the simulated stem and diffused passively throughout the stem. Its interaction with H resulted in the degradation of both H and M (Reaction X 11 : M+H!Ø).
Additional elaborations of the HBPM model. To further increase biological realism, a final iteration of the HBPM model added an additional molecule, A, which was produced exclusively at the top layer of the simulated stem (apex) and interacted with H in a regulatory feedback loop. In the augmented model, A, rather than H, was transported by P basipetally and upward facilitated transport did not occur. This augmented model did not involve any additional pattern formation mechanisms compared to the HBPM model, but it incorporated source-sink dynamics in plant vascular pattern formation in which the stem apex was the source of A. This augmented HBPM model is detailed in S2 Table. Examining hypotheses in silico with the HBPM model Effect of stem diameter on formation of vascular bundles. SVBs are most frequently found in wide parenchymatous stems, hypocotyls, or roots, and they tend to be absent from thin ones [29,30,47]. Can reductions in stem diameter-and in particular changes in size of the shoot apex near to where vascular patterning is initiated by procambia-explain the absence of SVBs in narrow stems and their development in wider stems? In narrow stems, M was predicted to diffuse into stem centers and suppress H (and therefore SVB) production throughout the stem. In contrast, H was predicted to be active in the centers of thick structures where M does not reach, since M is depleted (reactions X 10 and X 11 , Table 1) while it diffuses from the cortical site of production toward the stem center (Fig 2).
To investigate the effect of rates of diffusion and production of M on the localization of SVBs, two rates of diffusion and two rates of production were explored (S1 Table, Simulations 3-6, parameters D M and k 9 ). Simulations 3 and 7-9 (S1 Table) explored the impact of varying the stem radius while keeping all other rate and diffusion parameters constant.
Effect of diffusion rates of H and B on density and size of vascular bundles. To investigate how diffusion rates of H and B may influence size and density of SVBs, their diffusion rates were varied while other reaction rate parameters were kept constant. Following the results of Wei and Winter [37], it was expected that higher diffusion rates of D H would result  Table 1. Dark lines with circle end points indicate suppression in the direction of the circle. Arrowed lines indicate activation in the direction of the arrow. Arrows with only one connected end represent constitutive or background production (arrow pointing to box) or degradation (arrow pointing away from box). Squiggly arrows represent diffusion. Some simulations also included mutual inhibition between M and B, but this was removed in later simulations.
https://doi.org/10.1371/journal.pone.0219055.g004 3D model of stem vascular development in wider simulated vascular bundles and higher D B /D H ratios would tend to produce lower densities of simulated vascular bundles. For these simulations, the 2D HB model that did not include M and P was used (i.e., Eqs 1 and 2) in order to focus on the role that diffusion rates of H and B play in patterning. D B /D H ratios were set to 10, 100, 1000, 2000, 2500, 3000, 3500, and 4000, and D H was set to 0.016 or 0.05. The SSA was run three times for each combination of parameters, and the regions of high [H] were counted, measured, and graphed as a function of D B /D H vs. D H once stationarity was reached.
Simulations typically took several days to reach stationarity. 2D models run more quickly than 3D models and therefore permit a more rapid exploration of model parameter space. However, to check that the expectations concerning diffusion rate and size and density of simulated SVBs were met in a 3D model, three 3D simulations (Simulations 12-14, S1 Table)

HBPM modeling of stelar patterning
During the course of developing the HBPM model for the above simulations, multiple combinations of parameters were explored with the goal of determining the range of applicability of the model to the general patterning of vascular steles. Parameters used during these simulations (Simulations 15-18) are described in detail in S1 Table.

Model creation
In deterministic versions of the Gray-Scott-Schnakenberg model that produce stationary regions of high [H], these regions tend to be uniformly spaced. In the 3D HBPM model, P appeared to canalize and spatially fix the more randomly dispersed conduits of high [H]. Simulation 1 was run for 279,104,164,995 reaction events after which the pattern of [H] did not change substantially (Fig 3B-3D). In 2D simulations that lacked P, this canalization of high [H] did not occur, and regions of high [H] moved until uniformly spaced (e.g., S11 Movie). The original HBPM model was augmented to include A, representing auxin, which is apically produced and transported basipetally in the model. Results of two simulations of the augmented model are presented. S19 Movie presents the outcome of a simulation with apical production of A, basipital transport of A via P, and diffusion of H, whereas S1 Fig presents the outcome of the same augmented model, but with no diffusion of H.

Examining hypotheses in silico with the HBPM model
Effect of stem diameter on formation of vascular bundles. Simulations 3-6 explored the impact of M by varying its rates of production and diffusion within a constant diameter stem. These simulations all produced longitudinally oriented SVBs and zones that lacked SVBs that corresponded to zone B of Fig 2. The rate of production of M most highly impacted the width of zone B (Fig 3E-3L; S3-S6 Movies).
To examine the impact of changes in stem diameter on SVB patterning, reaction parameters were held constant in narrow, medium, and wide stems. Simulation of a wide diameter stem resulted in scattered SVBs that were restricted to the central part of the stem. There were a few, scattered SVBs in medium-diameter stems, and M suppressed the production of SVBs in narrow stems (Fig 3M-3O; S3 and S7-S9 Movies).

Effect of diffusion rates of H and B on density and size of vascular bundles.
Holding all other parameters constant while varying the diffusion rates of H and B resulted in changes to the number and widths of simulated vascular bundles. As predicted, changes to the diffusion rate of H (y-axis, Fig 5A) in large part accounted for variation in spot size, whereas variation in the ratio D B /D H (x-axis, Fig 5A) in large part accounted for changes in spot density (Fig 5; S10 and S11 Movies provide example 2D simulations with contrasting H and B diffusion parameter combinations). However, density and spot size were not independent, and spot sizes tended to decrease with increased spot density, even for constant diffusion rate of H and varying diffusion rate of B. Low diffusion rates of H resulted in spatially noisy [H] and no persistent regions of high [H] (S10 Movie).
All simulations examining diffusion rates of H and B were carried out in 2D did not include M and P. Without P, the longitudinal canalization effect was not present and spots spread uniformly throughout cross sections in the absence of M (S11 Movie). 3D simulations of the HDPM model with varying D H and D B matched expectations. Higher diffusion rates of H resulted in larger zones of high [H], and smaller spots tended to be in higher density than larger spots, even for larger D B /D H ratios. (S12-S14 Movies, S1 Table).

HBPM modeling of stelar patterning
The HBPM model recreated multiple stelar patterns with or without SVBs. A ring of vascular bundles (eustele) typical of the primary state of growth in eudicots resulted when B was degraded at a low background rate (Fig 3P; S15 Movie). When the degradation rate was further lowered and production of M decreased, a ring of radially-arranged vascular bundles with central SVBs formed (Fig 3Q; S16 Movie); the patterning of vascular bundles resembled the A (cortex), B (lack of SVBs) and C (presence of SVBs in pith) zonation of kohlrabi (Fig 2B). When background degradation of B was removed, a vascular cylinder (siphonostele) was produced with SVBs internally (Fig 3R; S17 Movie). The HBPM model established both the positioning of the secondary vascular cambium at the stem periphery at the interface between zones A and B, and it recreated the centrally-located SVBs. The zone B that lacks SVBs appeared to result from depletion of M and B in region B caused by the formation of the radially oriented bundles at the stem periphery. As a result, the core reaction that enhanced production of H (X 1 of Table 1) could not occur in zone B where [B] was depleted.
Additional stelar patterns were observed under different parameterizations. In the simulation of a narrow stem, the results resembled a central rod of vascular tissue (haplostele) ( Fig  3T; S8 Movie). Such a protostele is typical of roots with narrow root apices and represents another instance of how changes to boundary conditions influence vascular patterning. When suppression by M was reduced, a layout qualitatively similar to that of some monocots (atactostele) resulted in a uniform ring of high [H] that resembles vascular bundles peripherally, and scattered bundles internally (Fig 3V; S18 Movie). Finally, a deterministic version of the Gray-Scott model, Simulation B, produced patterns in simulated [H] that resembled plectosteles in stems of some lycopods (Fig 3W and 3X; S2 Model file for use with Ready 0.8 program).

Discussion
The HBPM model recreated three dimensional plant stem vascular patterning using a Turinglike mechanism in a stochastic, reaction-diffusion framework. Vascular patterns it recreated included atactostele, plectostele, haplostele, siphonostele (continuous vascular cylinder), and eustele (siphonostele with vascular bundles arranged radially) arrangements (Fig 3). Additionally, the model recreated patterning of a cambial variant in which supplemental vascular bundles (SVBs) develop in the center of the stem. As a mathematical model, it requires further grounding in molecular biological details, and it offers several testable predictions about evolutionary-developmental mechanisms of pattern formation.

Evolutionary-developmental hypotheses stemming from simulated outcomes of the HBPM model
Multiple sets of simulations explored how (1) changes in stem size (boundary conditions), (2) changes to the spatial expression of an inhibitory molecule (M), and (3) changes to diffusion rates of H (the molecule triggering vascular development) and B (its substrate) can influence vascular patterning.
(1) Changes to stem diameter altered SVB patterning in the HBPM model. Wide stems tended to develop longitudinally-oriented SVBs in their centers whereas narrow stems lacked SVBs. In these simulations, no changes to reaction, production, degradation, and diffusion rates were made; only stem diameters were altered. Empirically, Torrey [48,49] found that stelar patterns in pea roots are largely a function of the diameters of apical regions.
One implication of the simulations that modeled changes to stem diameter is that evolutionary gains or losses of a trait may have little to do with evolutionary changes to gene expression rates or functions of the genes responsible for the development of a trait. Instead, losses, for example, may be due to evolutionary changes in boundary conditions that allow a diffusible inhibitor to reach the location of pattern formation that was previously inaccessible to the inhibitor. This scenario can explain the loss or gain of SVBs in the HBPM model when stem diameter is varied. When (simulated) evolution decreases stem girth (through other developmental pathways), SVBs fail to develop due to inhibition of H by M, whereas (simulated) evolution of wide-stemmed plants in which the inhibitor molecule does not reach the centers of stems will result in SVB development (compare Fig 3M-3O). So, changes to genes responsible for the evolutionary loss or gain of a trait may be independent from the molecular-developmental pathways responsible for the trait's development.
(2) Changes in spatial expression of the inhibitor, M, altered the position of the stele and vascular bundles in relation to the edge of the stem (Fig 3E-3L). With no expression of M, the atactostele pattern found in wide-stemmed monocots was the result. Evolution can cause such changes by altering promotor regions, thereby influencing the rate and spatial positioning of expression. Carroll and others [1,50] argue that many changes to organismal patterning are due to such mechanisms.
(3) Changes in diffusion rate of H and B altered the number and size of SVBs even when all other aspects of the model (e.g., reaction kinetics, production and degradation rates, etc.) remained constant (Fig 5). Moreover, simulations in which diffusion rates of H fell below a threshold rate failed to develop vascular bundles altogether (S10 Movie). Mathematical theory suggests that inhibitor molecules must diffuse at least six to seven times faster than activator molecules in order for heterogeneous patterns of the activator molecule to emerge ( [51] and p. 246 in [52]).
A broader hypothesis that follows from this modeled result is that phenotypic diversification may occur as a consequence of evolutionary changes to protein movement rates. The logic of this hypothesis is as follows: (1) Evolution can influence lengths of proteins through insertion or deletion (indel) mutations; (2) changes in protein lengths influence their rates of movement (smaller molecules tend to move faster); (3) changes in rates of molecular movement can impact patterning processes during development. Therefore, evolution may be able to alter patterning by changing molecular movement rates. The first statement follows by definition of indel mutations. The second statement follows from basic thermodynamics and is empirically supported by some relatively recent work that tracked individual, tagged molecules. This research found that protein size is the primary determinant of cytosolic diffusion rate in E. coli cells [53]; Nenninger et al. [53] also found that, empirically, the coefficient of diffusion is related to the size of the protein under constant temperature and fluid viscosity via the Stokes-Einstein formula. Finally, the third statement has been recognized since Turing's 1952 paper [25]. If molecular movement rates are largely a function of molecular size, then even in the absence of selection for change in amino acid composition and protein function, variation in protein size can play a role in phenotypic evolution through changes in movement dynamics. This relationship between protein size and patterning of serially homologous features also raises the intriguing possibility that different alternative splice isoforms of a gene's transcript permit fine tuning of diffusion rates-and hence pattern formation processes-in specific cell lineages.
Interestingly, Lipman et al. [54] discovered that length variation is highest among the most conserved, "more important" (their words) proteins. They interpreted protein size variation in terms of transcriptional costs, and they also suggested that size variation in conserved proteins reflects divergence in protein function. Selection on protein diffusion rate provides an alternative interpretation in which the fine tuning of protein length (hence diffusion rate) and not protein function, per se, influences the outcomes of spatio-temporal developmental patterning and selection on phenotypes.

Biological assumptions and possible molecular bases of the HBPM model
In order to design empirical experiments to test the above hypotheses based on the HBPM model, the components of the model must be empirically grounded with molecules of relevance to vascular patterning in plants. Mathematical modeling and molecular-genetic research about plant vascular development have long suggested that members of the Class III HOMEO-DOMAIN-LEUCINE ZIPPER (HD-ZIP III) gene family are essential for the initial specification of growth regions (procambia) that develop into xylem and phloem [5,[55][56][57][58][59][60][61] through regulatory feedback with auxin [62], KANADI [60], brassinosteroid (BR) [4,63], and possibly with Class I KNOTTED-LIKE HOMEOBOX (KNOX) transcription factors [64].
The model of Benítez and Hejátko [5] required positive autoregulation of an HD-ZIP III gene (HB8 in Benítez and Hejátko) in coordination with another molecule for vascular patterning to emerge. A similar model of autocatalysis has been proposed for at least one other HD-ZIP III gene family member, PHABULOSA (PHB), in which a ligand for PHB activates the PHB protein and may increase its expression [65]. Empirically, the HD-ZIP III gene HB8 is upregulated in regions of tuberous stems and roots with SVBs in Brassica crops [36]. In the HBPM model, H, appears to play the role of an HD-ZIP III protein involved in procambial specification and xylem differentiation, whereas B plays the role of the ligand involved in H's autoregulation (Table 1 Reaction X 1 ). HD-ZIP III transcription factors are active as homodimers [66] consistent with H's dimeric interaction in the core HB model (Eq's 1 and 2). B remains an unspecified positive regulator of the HD-ZIP III protein that is itself degraded by the HD-ZIP III protein. This regulator of the HD-ZIP III gene may be BR, an ARF such as MP, auxin, a KNOX I transcription factor, or an unidentified molecule [67]. Ibañes et al. [4] suggested that BR is essential for vascular bundle formation in Arabidopsis, and Benítez and Hejátko [5] postulated that HB8 self-upregulates via BR. In the current model, B may play the role of BR or a regulator of BR such as a member of the BRL gene family [63].
Ibañes et al. [4] also suggested that polar auxin transport is required for the patterning of vascular bundles. Sachs [68][69][70] was the first to suggest that the self-reinforcement of auxin through polar transport canalizes auxin at the site of procambial strand initiation. Basal positioning of the auxin efflux protein PIN1 within cells via gravity-sensing mechanisms [71][72][73] precedes this canalization effect [21,60,74,75]. Specification of procambial cells occurs through the activities of auxin response factors (ARFs) such as MP and the induction of HD-ZIP III HB8 expression [62] in a feedback loop. Moreover, in the model of Muraro et al. [7], an HD-ZIP III gene (PHABULOSA) was an indirect positive regulator of PIN proteins, and ARFs were indirect positive regulators of HB8 in Benítez and Hejátko's [5] model. Furthermore, HD-ZIP III genes can promote expression of the PIN-FORMED (PIN) efflux carrier proteins [62,74] that contribute to polar auxin transport in a complex network involving positive regulation of HD-ZIP III transcription factor genes (e.g., HB8) by the auxin response transcription factor (ARF), MONOPTEROS (MP) (reviewed by, e.g., [12]). During simulations of the HBPM model, no parameter combinations were found that recreated the patterning of longitudinally-oriented vascular bundles in 3D without simulating an efflux molecule, P. Based on the above considerations, P matches the function of PIN1.
Lastly, how does M function in the HBPM model? Empirically, microRNA 165/166 depletes HD-ZIP III and thereby limits the development of vascular tissue [76]. M performs the function of microRNA 165/166 in the HBPM model. In roots, SHORT ROOT (SHR) and SCARE-CROW (SCR) that accumulate in the endodermis increase microRNA 165/166 production in the endodermis (reviewed by [77]). In the HBPM model, expression of M is restricted to the stem cortex a priori, analogous to how its expression is restricted to the endodermis in roots [7]. However, the empirical spatial distribution of microRNA 165/166 in stems is less well understood. The model suggests that miRNA165/166 suppression of an HD-ZIP III protein may restrict the differentiation of SVBs to the central regions of plant structures and foster the partitioning of stems into the familiar cortex-cambium-xylem-pith configuration.

Addressing biological assumptions
If, indeed, H, B, P, and M correspond to an HD-ZIP III protein, a brassinosteroid (or unknown molecule), a PIN efflux protein, and microRNA 165/166, respectively, then the initial HBPM model needs refinement to better correspond to biological reality. First, PIN1 transports auxin and is very unlikely to transport HD-ZIP III proteins. Second, auxin is known to be produced apically and it is transported primarily basipetally. Third, HD-ZIP III proteins are quite large (~91kDA), and it is unlikely that they can diffuse freely between cells, even through plasmodesmata [78]. To increase biological realism, HBPM was augmented to model auxin, symbolized as A, explicitly. The model was reconfigured so that A, and not H, was transported basipetally by P. The simulated expression of A was also restricted to the top layer of the simulated stem to model apical expression of auxin. H and A were simulated to regulate each other in a feedback loop in the augmented model. These additional reactions are presented in S2 Table. In the augmented model simulations, the concentration of [H] remained high in the apical region. Below this region of high concentration was a region of low concentration, and under this region, longitudinally oriented regions of high [H] occurred matching the layout of procambial patterning in some plants (the discontinuity in S19 Movie corresponds to a shift of view from the top layer of high concentration to the lower layer in which longitudinal patterning of high [H] developed).
When the diffusion rate of H was set to 0 to limit its diffusion due to its large size, longitudinally arranged regions of high [H] resulted, but these were irregularly spaced in the stem and had narrow diameters. In all augmented simulations with A, longitudinally arranged regions of high H were crooked, unlike those in the "un-augmented" HBPM simulations.
So, the HBPM model captures a lot of the variation and patterning in plant stem vascular tissue, but it remains incomplete. The actual regulatory network underlying vascular development is substantially more complicated than a four-or five-molecule model, and additional components are required to further refine the model of three-dimensional vascular patterning.
In particular, longitudinal, cellular growth near the apex and in rib meristems was ignored. Gene expression states of cells could be epigenetically inherited during cell division in these regions resulting in longitudinally oriented vessels without PIN proteins. Future modelling efforts should include growth dynamics to examine how apical growth may contribute to three dimensional patterning of vascular tissue in stems.

HBPM model predictions
Further functional molecular genetic validation of the HBPM model is required. This model makes multiple predictions that can be tested in the wet-lab: 1. Production of miRNA165/166 is highest at the stem periphery, and production of miRNA165/166 may be reduced in monocot stems with atactosteles.
2. Release of HB8 or other HD-ZIP III genes from suppression by miRNA165/166 will increase the zone of SVB production throughout the stem in stems composed of ground tissue or parenchyma, more generally (c.f. [58]).
3. Suppression of HB8 expression (or possibly a different HD-ZIP III paralog) will annihilate SVB development.
Different parameterizations of the HBPM model led to patterning of vascular bundles that resembled monocot vascular bundle patterning (Fig 3V) or eudicot patterning ( Fig 3P) and recreated the partitioning of stem into cortex, vascular cambium, and pith. M influenced the relative sizes of these regions, so altering the location and expression rate of microRNA 165/ 166 is also predicted to influence the relative sizes of cortex, vascular regions, and pith in the stem.  Table 1) accounts for much variation in plant stem vascular structure. In each panel of the figure, the top represents the shoot apex, followed by a cross section with the primary state of growth (compare to Fig 3P), followed by a cross section of the mature stem. (A) Monocot stem vascular patterning as seen in corn or bamboo (compare to Fig 3V). (B) Woody plant with primary vascular bundles that merge and develop into the vascular cambium which produces lignified xylem (black sectors). Reduced axial parenchyma and pith may impose limits the production of SVBs although these are occasionally found in pith of woody plants [27]. (C) An herbaceous eudicot maintains the primary state of growth with limited cambial development and limited increase in stem girth. Narrow stems restrict the production of SVBs (compare to Fig 3O and 3P). (D) A stem-succulent plant with extensive parenchymatous wood. The primary vasculature transitions into a secondary vascular cambium that produces parenchyma primarily. H and B interact in wide zones of parenchyma to produce SVBs in the stem centers away from the suppressive effects of M (compare to Fig 2B and Fig 4Q-4S).

Conclusions
The HBMP model and its augmentations demonstrate that, in principal, an activator-substrate-inhibitor Turing-like patterning mechanism with longitudinal transport via Michaelis-Menten kinetics of a regulatory molecule is sufficient to canalize the differentiation of radiallypatterned, longitudinally-oriented vascular bundles within stem ground tissues. The result that a single reaction network can recreate most of the common vascular patterns in plant stems raises the possibility that a single conserved genetic network regulates stem vascular patterning in plants more broadly and that evolution may "configure" this network within particular lineages by altering boundary conditions and reaction-diffusion parameters.
Although incomplete, the augmented HBMPA model appears to correspond best with activities of HD-ZIP III regulatory gene products (H), miRNA inhibitors of H (M), PIN1 efflux carrier proteins (P), auxin (A), and an unknown substrate (or cofactor) of H (hypothesized to be a brassinosteroid, B). The HD-ZIP III genes comprise an ancient family of genes that are active in nonvascular and vascular land plants, and they appear to be involved in vascular development broadly throughout the vascular plants. In animals, it is increasingly becoming apparent that evolutionary tweaking of the temporal and spatial expression of conserved "toolkit" genes [50] underlies animal diversity. Is there a "toolkit" for vascular development in plants? If so, HD-ZIP III gene products and molecules interacting with them are likely to be important components of the core "kit".
The current model provides a unified framework that suggests that shared, homologous, deeply conserved molecular genetic machinery is sufficient to explain diverse stem vascular patterns (Fig 6). However, the molecular genetic details of the HBPM model and its evolutionary-developmental implications concerning the fine-tuning of biological patterning via evolution of diffusion rates and reaction-diffusion boundary conditions require further empirical validation and scrutiny. (VTI) S1 Program. Executable Java archive file (.jar) to run simulations. This file is a ZIP archive file of multiple files. In addition to the executable jar file, the ZIP file includes a tutorial, source code in the .jar file, and java library files (lib directory) required to run the reaction-diffusion program. (ZIP) S1 Table. Simulation parameters. Column headers indicate the simulation lattice dimensions, the reactions (from Table 1) and their rate parameters, and the diffusion parameters. S1-S18 Movies file names and figure numbers are also provided to see visualizations of the simulations. The columns "radius M" and "width M" describe the radius and thickness of the cylindrical shell in which M is constitutively expressed. (XLSX) S2 Table. Augmented HBPM model with apical production of A and model parameters. This table includes parameters for S19 Movie and associated simulation. (XLSX) S1 Movie. Simulation 1 Three dimensional simulation of the HB model. Movie files are Quick-Time-formatted movie files that visualize the time course of simulation data. Parameter values for corresponding simulations are in S1 and S2 Tables. The figure number that visualizes the end result of the simulation is provided next to each movie file name. In all movies, the top, middle, and bottom row of each column represent the concentrations of reactants in cross-section, longitudinal section, and total concentration over time, respectively. Color bars follow conventions described in