Gillespie-Lindenmayer systems for stochastic simulation of morphogenesis

Lindenmayer systems (L-systems) provide a useful framework for modelling the development of multicellular structures and organisms. The parametric extension of L-systems allows for incorporating molecular-level processes into the models. Until now, the dynamics of these processes has been expressed using differential equations, implying continuously valued concentrations of the substances involved. This assumption is not satisfied, however, when the numbers of molecules are small. A further extension that accounts for the stochastic effects arising in this case is thus needed. We integrate L-systems and the Gillespie’s Stochastic Simulation Algorithm to simulate stochastic processes in fixed and developing linear structures. We illustrate the resulting formalism with stochastic implementations of diffusion-decay, reaction-diffusion and auxin-transport-driven morphogenetic processes. Our method and software can be used to simulate molecular and higher-level spatially explicit stochastic processes in static and developing structures, and study their behaviour in the presence of stochastic perturbations.


Introduction
L-systems are a powerful formalism for modelling the development of growing linear and branching structures, from basal filamentous organisms to trees and entire plant ecosystems (Prusinkiewicz and Lindenmayer 1990;Lane and Prusinkiewicz 2002). The introduction of parameters into L-systems (Lindenmayer 1974) has made it possible to capture genetic regulatory mechanisms and physiological processes that affect development by using a continuous-level representation of the substances involved (Coen et al. 2004;Allen et al. 2005;Buck-Sorlin et al. 2005;Prusinkiewicz et al. 2009). However, in many biological processes the numbers of molecules are relatively low, and the assumption of continuity of concentration is not satisfied (McAdams and Arkin 1997;Elowitz et al. 2002;Rao et al. 2002). In these cases individual molecules must be considered, and algorithms that simulate stochastic fluctuations in molecule numbers are needed.
Stochastic simulation of reaction kinetics was pioneered by Gillespie (1976Gillespie ( , 1977. The essence of his method is a discrete-event simulation of reactions between individual molecules in a spatially homogeneous, well-mixed system. To simulate heterogeneous systems, Gillespie (1976) proposed an extension, in which molecules diffuse between subvolumes in compartmentalized space. Spicher et al. (2008) incorporated Gillespie's algorithm into P-systems (Păun and Rozenberg 2002), where the modelled system, such as a cell, is structured into a hierarchy of compartments enclosed by 1. determine the delay τ (inter-reaction time) with which the next reaction will take place, and the index of the next reaction j ∈ (1, . . . , M), 2. modify the state X, taking into account the reactants removed from the system and products added to the system by reaction R j , and 3. advance simulation time t by τ .
To generate one random pair (τ , j), Gillespie proposed two methods: the first-reaction method and the direct method. In the first reaction method, a putative time τ j , j ∈ (1, . . . , M), is generated for each reaction, and the reaction with the smallest time, τ j , is chosen. In the direct method values for τ and j are generated according to a joint probability function of τ and j. The direct method is more efficient than the first reaction method because fewer random numbers must be generated per iteration step. Consequently, we only consider the direct method.
Gillespie showed that the time between two reactions can be described by an exponential distribution. First, the combined propensity of all possible reactions is computed: (1) The inter-reaction time τ is then calculated as an exponentially distributed random variable with a mean value of 1/α 0 (X). The index of the next reaction is described by a discrete probability distribution, where α j (X)/α 0 (X) is the probability that the next reaction is R j . The interreaction time τ and the next reaction R j are then generated according to these probability distributions using the inversion method (Gillespie 1976). Specifically, given two independent random numbers r 1 and r 2 , generated with uniform distribution in the interval (0, 1], the interreaction time is obtained using the formula Table 1. Terms involved in the calculation of reaction propensities for several reaction types. Column h: the number of distinct molecular combinations of reactants S 1 and S 2 as a function of the number or their molecules, X 1 and X 2 . Column c: the stochastic reaction parameter as a function of the reaction rate k. Symbol ∅ denotes an external source of products. Parameter Ω is the volume in which the reactions take place.

Reaction name Reaction h c
Source ∅ k → products 1 kΩ and the reaction index j is determined by solving the equation To incorporate spatial inhomogeneities into simulations of chemical systems, Gillespie (1976) proposed an extension to the basic Stochastic Simulation Algorithm, in which the volume V is divided into n subvolumes (also referred to as compartments or components) V i (i = 1, 2, . . . , n). Two basic ideas underly this extension: • Propensities of reactions taking place in different subvolumes are calculated individually for each subvolume. Thus, instead of a reaction propensities α j , the algorithm calculates propensity α ij of reaction j taking place in subvolume i. • Transport of a molecule from subvolume i to its neighbour î is treated as a unimolecular reaction that removes a molecule from subvolume i and deposits it in subvolume î.
The inter-reaction time τ , the index j of the next reaction or transport event, and the index i of the subvolume where the next event occurs-or subvolume pair (i, î ) for transport events-are then generated using Equations (1-3), in which the summation extends to all pairs (i, j) of reaction events and triplets (i, î, j ) of transport events. Formal details are presented by Stundzia and Lumsden (1996).
With the extension to subvolumes, Gillespie's algorithm can simulate the time evolution of a spatially explicit system in which the partition into subvolumes is fixed, while the molecules within each subvolume are distributed uniformly. To model the development of multicellular organisms, however, it is necessary to consider spatial structures in which the number of subvolumes and reactions associated with them may change over time (e.g. following cell division). We integrate Gillespie's algorithm with L-systems to provide a mechanism in which both fixed and developing structures can be simulated easily.

L-systems
Parallel rewriting systems, subsequently called L-systems, were introduced by Lindenmayer (1968) to specify, model and reason about the development of multicellular structures (whole organisms or their parts) with filamentous or branching topology. A structure is represented by a string of symbols (letters) that correspond to its individual components, such as cells, higher level architectural units or compartments resulting from a discretization of a continuous space. The evolution of the structure state over time is characterized by rewriting rules, also called productions, which specify how a predecessor symbol is replaced by zero, one or more successor symbol(s) in the string. For example, the rule A → BC may be used to represent the division of cell a into two daughter cells, B and c. The rules are applied in parallel to the entire string, to capture the simultaneous progress of time in all parts of the organism (Prusinkiewicz and Lindenmayer 1990).
Parametric L-systems associate numerical parameters with the symbols (Lindenmayer 1974;Prusinkiewicz and Lindenmayer 1990;Prusinkiewicz et al. 2018), for instance to quantify the chemical substances in each component of the structure. For example, the rule C(z) → C(z − µ z z∆t) describes the decay of substance Z in cell c, where z is the concentration of Z at time t, µ z is the decay rate and ∆t is the time step.
A symbol with the associated parameters is referred to as a module. Communication between modules and the transport of substances within the structure can be modelled using context-sensitive productions, in which the production successor depends not only on the predecessor module, but also on its neighbours or context. Notationally, the context is separated from the predecessor by symbols < and >. For example, production specifies that concentration a of some substance in cell c changes according to the fluxes J l and J R through walls W on both sides of c. Context-sensitive productions facilitate modelling of multicellular structures, because they eliminate the need to index cells, and then reindex them as the structure develops and neighbourhood relations change.
We will specify L-systems using the L+C modelling language implemented in the simulation program lpfg, which augments the expressive power of L-systems by combining them with C++ (Karwowski 2002;Karwowski and Prusinkiewicz 2003;Prusinkiewicz et al. 2007Prusinkiewicz et al. , 2018. L+C is relatively self-explanatory; for example, the context-sensitive production (4) is written in L+C as In a complete L+C program, symbols representing components of the structure are declared using the keyword module, with the parameter types listed in parentheses. The initial string is specified after the keyword axiom. Predefined keywords Start and End indicate optional blocks of C++ statements executed at the beginning or end of the simulation, for example to initialize variables used in the simulation or report statistics gathered during its execution. Likewise, StartEach, and EndEach indicate statements executed at the beginning or end Downloaded from https://academic.oup.com/insilicoplants/article-abstract/1/1/diz009/5614990 by University of Calgary user on 02 December 2019 in silico Plants https://academic.oup.com/insilicoplants of each simulation step. Following these blocks, the L-system productions are listed. A production may have more than one successor, each preceded by the produce keyword. In the original version of L+C the applicable successor is selected using a conditional statement; in Gillespie L-systems the selection mechanism is extended to include the stochastic mechanism described in the next section. The set of all productions can be partitioned into subsets called groups by the statement group: id (Prusinkiewicz et al. 2007), interspersed between productions. This statement assigns a numerical identifier id to all productions following it in the production list, until the next group: id statement occurs or the production list ends. The group: id statements are used in conjunction with the function UseGroup(id), which is typically called within the StartEach block and specifies which group of productions should be used in the forthcoming simulation step. The notion of groups thus increases the flexibility of L-system programming by making it possible to employ different production groups in different simulation steps.

Integration of Gillespie's algorithm and L-systems
To combine Gillespie's algorithm and L-systems, we introduced the notion of a Gillespie group of productions, identified by the keyword ggroup. A Gillespie L-system may include both 'ordinary' and Gillespie groups, but only one group is active in any simulation step. Productions in an ordinary group are applied to all modules in the string in parallel, consistent with the standard definition of L-systems. In contrast, a single module in the string and a single production or production successor are selected in a simulation step using a Gillespie group. This selection is effected using Gillespie's Stochastic Simulation Algorithm extended to subvolumes-identified with modules-with probabilities controlled by the expressions specified after the propensity keyword associated with each successor. As the selected production is applied, the current simulation time is advanced by the inter-reaction interval. The current time (integrated from the beginning of the simulation) can be accessed using the predefined L+C function GillespieTime(). Time management and event scheduling may require attention in simulations combining stochastic and deterministic productions; an example is discussed in Auxin-driven morphogenesis section (Program 4b).
The following example illustrates the method. Consider substances A and B decaying in two types of modules, C and D, declared as follows: The decay events are defined by productions in the Gillespie group (although it is the only group, it must be selected explicitly using the UseGroup(1) statement): In general, the propensity of each reaction is the product of a stochastic reaction parameter and the number of combinations of the reacting molecules. As decay is a unimolecular reaction, the stochastic reaction parameter is equal to the reaction (decay) rate µ, and the numbers of combinations are equal to the numbers of reactant molecules (Table 1). Suppose the rates are µ a = 0.1 and µ b = 0.2, and the current state of the system is described by the string of four modules A simulation step begins with the simulator (lpfg) calculating the propensity of each applicable reaction in each individual module, producing the results shown in Table 2. On this basis, the simulator stochastically selects a single reaction taking place in a single module using Equation (3). For instance, the probability of selecting reaction A → ∅ in module C is 0.4 2.0 = 0.2, the probability of selecting reaction B → ∅ in module C is also 0.4 2.0 = 0.2 and the probability of selecting reaction A → ∅ in the last module D is 0.5 2.0 = 0.25. Once the reaction and module have been selected, the simulator decrements the number of molecules of the reacting substance in the affected module by one, as specified by the production successor, and determines stochastically the inter-reaction time using Equation (2) with α 0 = 2.0. The next simulation step can then be performed, beginning with the recalculation of propensities. The entire simulation run is a sequence of such steps.

Petri nets
As the reactions become more complicated, it is convenient to represent them graphically using stochastic Petri nets (Goss and Peccoud 1998). We use them in this Table 2. Example of propensity calculations for two decay reactions in a string of four modules. paper to visually complement textual specifications of reaction systems. A Petri net is a directed graph with two types of nodes: places, drawn as circles, and transitions, drawn as rectangles (Fig. 1). In the context of chemical processes, places represent substances, and transitions represent reactions or movement of molecules between components. The nodes are connected by directed arcs (arrows), such that the arcs pointing from places to transitions denote the reactants, and arcs pointing from transitions to places denote the products. Arcs connecting places to places or transitions to transitions are not allowed. In a stochastic Petri net transitions are executed or fired by a stochastic process: in the scope of this paper, using Gillespie's Stochastic Simulation Algorithm.

Examples
We now illustrate Gillespie L-systems with a sequence of examples of increased complexity.

The Lotka-Volterra process
The first example revisits the Stochastic Simulation Algorithm implementation of Lotka's (Lotka 1920) chemical process with an oscillatory behaviour, originally presented by Gillespie (1977). The purpose of this example is to show the basic structure of an L+C program using a Gillespie group of productions. The reactions are: where Y 1 and Y 2 are two chemical substances, and k 1 , k 2 and k 3 are the reaction rates. The stochastic Petri net describing this system is shown in Fig. 2, and the corresponding L+C implementation is given in Program 1.
Program 1. L+C implementation of Lotka's chemical system.
The program operates on a single module C representing the entire reaction volume, with two parameters specifying the current number of molecules Y 1 and Y 2 . We assume that the volume Ω in which the reactions take place is equal to 1, and thus the stochastic reaction parameters c j for all reactions, including the bimolecular reaction R 2 , are equal to their reaction rates k j (Table 1). In each simulation step, the simulator computes the propensities of the applicable productions (in this example, all three of them are always applicable), selects one using Gillespie's direct method, applies it to module C and advances the current simulation time by the inter-reaction interval. A sample run of Program 1 is shown in Fig. 3. Consistent with the theory (Lotka 1920), the number of molecules exhibits an oscillatory behaviour. The jagged character of the plots, clearly seen in Fig. 3B, reflects the inherently discrete and stochastic nature of the Lotka system, as the molecules only occur in integer numbers and enter into reactions randomly.
The labels a 0 , B 0 , c 0 and D 0 indicate the initial (t = 0) numbers of molecules of each substance. Independently of Lotka's chemical system, Volterra (1928) proposed an ecological model of fish catches in the Adriatic Sea that exhibits the same oscillatory behaviour. In this case, the reactions are interpreted as follows: R 1 : a prey species Y 1 reproduces while feeding on some food source that does not deplete over time, R 2 : a predator species Y 2 reproduces while feeding on the prey species, and R 3 : the predator species dies by natural causes.
This reinterpretation shows that applications of Gillespie L-systems are not limited to molecular-level simulations.

Diffusion and decay
Let us now apply a Gillespie L-system to model a simple spatially explicit process, in which substance A diffuses and decays in a one-dimensional medium. The standard description of this process has the form of the partial differential equation where a = a(x, t) is the concentration of A at point x and time t, µ a is its decay rate and D a is the diffusion rate (Edelstein-Keshet 1988). This model can be spatially discretized into a linear structure (a one-dimensional cell complex) with two types of components: cells and cell walls (Prusinkiewicz and Lane 2013). Following the law of mass conservation, the rate of change in concentration a i of substance A in cell i is then equal to: where J (i−1)→i , i = 2, 3, . . . , n − 1, is the flux of A through the wall between cells i − 1 and i. According to Fick's law, this flux is proportional to the concentration difference a i−1 − a i : A Petri net corresponding to Equations (6) and (7) is shown in Fig. 4. It leads to the Gillespie L-system implementation in Program 2a. The program begins with the definition of parameters: diffusion rate D a , decay rate µ a and the number X of molecules A in the boundary cells. Declarations of two module types-cell C and wall W-follow. The parameter of cell C is a non-negative integer denoting the number of molecules of substance A in this module. The parameter of wall W is the integer −1, +1 or 0, indicating whether a molecule will be transported through W to the left, to the right or not at all. Program execution is controlled by the StartEach and EndEach statement blocks, which alternate between ggroup 1 and group 2 in consecutive simulation steps. The Gillespie group, ggroup 1, has two rules. The first rule, with the predecessor W(dir l ) < C(a) > W(dirr), describes the decay of a single molecule of A. The context is included to maintain boundary conditions: the first and the last cell lack an incident wall (see the axiom), which fixes the number of molecules of A in them. The second rule, with the predecessor C(al) < W(dir) > C(ar), specifies diffusion of a molecule to the left or to the right using alternative propensity ... produce statements. The actual transport is effected by the first production in the standard L+C group 2, which operates in parallel on all modules C in the string except for the boundary cells. If the transport direction is left to right (dir l or dir r is 1), the cell to the left of the wall will lose a molecule and the cell to the right will gain one. Conversely, if the direction is right to left (dir l or dir r is −1), the cell to the left of the wall will gain a molecule and the cell to the right will lose it. The last production resets all diffusion events to 0, in preparation for the next iteration of the simulation.
Program 2a. Gillespie L-system implementation of the diffusion-decay process using explicit representation of walls. Program 2b. Alternative Gillespie L-system implementation of the diffusion-decay process using productions operating on cell pairs. Fig. 5A-C. The initial concentration of molecules A in the interior cells was assumed to be 0. We varied the number X of molecules of A in the boundary cells between different simulation runs to show the effects of increasing this number on the results. In each simulation run, the cell volume Ω was set to be numerically equal to X, so that concentration X Ω was equal to 1. This normalization facilitated comparisons of simulations with different molecule numbers. As expected (Gillespie 2007;Wang et al. 2007;Vigelius and Meyer 2012), the stochastic solution became less noisy as the number of molecules increased. For a comparison, Fig. 5D shows a solution to a deterministic diffusion-decay system with continuous representation of concentrations (Equations (6) and (7)). The apparent convergence to the deterministic solution is consistent with the convergence of the molecular motion produced by voxel-hopping to the standard diffusion equation (Gillespie et al. 2014). The stochastic simulation has the advantage of better representing the diffusion-decay process when the number of molecules is small.

Results of sample simulation runs are shown in
Incidentally, the same diffusion-decay process can also be simulated using the Gillespie L-system given by Program 2b. In this case, diffusion is effected by a single production (with two alternative produce statements) operating on a pair of modules. The boundary conditions, maintaining constant number of molecules in the first and last cell, are enforced by representing these cells using distinct modules B. The implementation of diffusion in a single simulation step makes Program 2b somewhat faster than Program 2a (~20 % in our implementation). On the other hand, Program 2a emphasizes the local character of diffusive transport, and is consistent with the standard definition of L-systems, according to which each production has a single predecessor (this assumption is key to the parallel operation of standard L-systems). The choice between Programs 2a and 2b is thus largely a matter of programming style.  (6) and (7)

Reaction-diffusion
In this example, we construct a Gillespie L-system to simulate a stochastic reaction-diffusion patterning process (Turing 1952;Gierer and Meinhardt 1972;Meinhardt 1982).We focus on pigmentation patterning in seashells, for which models expressed using partial differential equations are well understood Klingler 1987, 1988;Fowler et al. 1992;Meinhardt 2009). The activator-substrate variant of these models is described by the following equations (Meinhardt 2009): Each of these equations extends the diffusion-decay system, discussed previously, with terms representing reactions between activator A with concentration a and substrate S with concentration s. At the molecular level, these reactions have the form: Following Table 1, their propensities are α 1 = σΩ, α 2 = ρρ 0 S and α 3 = ρ Ω 2 SA(A − 1), respectively. The resulting stochastic Petri net combines these reactions with two diffusion-decay models: one for activator A and another for substrate S (Fig. 6). Correspondingly, Program 3, specifying the activator-substrate process in L+C, has production groups similar to Program 2. Figure 7A-C show the results of three runs of the simulation of the pigmentation pattern found in the seashell amoria undulata (Fowler et al. 1992;Meinhardt 2009). The images represent consecutive states of the simulation in a row of n = 100 modules, obtained for three different values of (mathematical) cell volume Ω: 10, 100 and 1000. As in the diffusion-decay model, the volume Ω was numerically equal to the initial number of molecules A and S in each run. The solutions used the following parameter values: ρ = 0.1 , ρ 0 = 0.005, µ a = 0.08, D a = 0.004, µ s = 0 and D s = 0. The σ parameter was modulated for each cell according to a sine function in order to generate lines of undulating shape (Fowler et al. 1992). Specifically, σ = σ min + (sin(2π · 3i/n) + 1)(σ max − σ min )/2, with σ min = 0.02 and σ max = 0.032, for each cell i = 1, . . . , n . For a comparison, Fig. 7D shows a numerical solution of Equations (8) and (9) assuming a continuous representation of concentrations. In nature, these temporal progressions take place on the growing shell margin, leaving a pigmentation pattern 'frozen' on the shell surface. The patterns found in a. undulata seashells (e.g. Fowler et al. 1992, Fig. 12) exhibit irregularities that are consistent with the stochastic simulation in Fig. 7C. The stochastic model thus gives a satisfactory explanation for irregularities that are observed in the natural patterns, but are not captured by the deterministic model. Program 3. Gillespie L-system productions implementing a stochastic model of the activator-substrate process.

Auxin-driven morphogenesis
In this last example, we consider the regulation of leaf shape by a key component of plant morphogenesis, the hormone auxin (Sachs 1991;Zažímalová et al. 2014). The export of auxin from a cell relies on the activity of carriers in the cell membrane. Among them, the PIN1 protein appears to play the most prominent morphogenetic role. The allocation of PIN1 to different regions of the membrane is regulated by auxin itself (Paciorek et al. 2005), creating a feedback loop that plays an essential role in many aspects of plant morphogenesis. Molecular-level details of this process are the subject of ongoing research (Abley et al. 2013;Cieslak et al. 2015), but within the leaf margin the end result is a preferential allocation of PIN1 to regions of the cell membrane abutting neighbouring cells with a high concentration of auxin (Hay et al. 2006;Scarpella et al. 2006). As a result of this 'up-the-gradient' allocation (Jönsson et al. 2006;Smith et al. 2006), a pattern of auxin concentration maxima and minima emerges. The maxima promote an outgrowth of future serrations, lobes or entire leaflets, thus shaping the developing leaf (Hay et al. 2006;Bilsborough et al. 2011, Bar andOri 2014;Runions et al. 2017;Conklin et al. 2019). The first model of the above process was formulated in terms of differential equations (Bilsborough et al. 2011).
Here we construct a stochastic model paralleling the version described by Prusinkiewicz and Lane (2013). The leaf margin is represented as sequence of cells that grow and divide upon reaching a threshold length. The volume of each cell is computed dynamically as the product of its cross-sectional area s, assumed to be constant, and length x, affected by growth and divisions. Within each cell, the model accounts separately for the number of PINs in the cytoplasm (PIN) and in the membrane regions abutting the left and right neighbouring cells (PIN l and PIN r , respectively). We distinguish between molecule count, denoted without brackets, for instance a for auxin molecules A, volumetric concentration [A] = A xS , and-in the case of molecules allocated to membrane regions-area concentrations, for instance [[PIN l ]] = PIN l S . PIN concentration on the membrane is the result of two processes: exocytosis, or the allocation of PIN from the cytoplasm to the membrane, and endocytosis, or the return to the cytoplasm. Consistent with the up-the-gradient polarization model, we assume that the propensity of exocytosis PIN → PIN l , which allocates a PIN molecule in cell i to the membrane region abutting its left neighbour i − 1, is proportional to the region area s, the PIN concentration in the cytoplasm, [PIN], and the auxin concentration in the neighbouring cell, [a l ]: We further assume that endocytosis PIN l → PIN, which deallocates a PIN molecule from the membrane to the cytoplasm, has a propensity characteristic of a decay process: Analogous equations apply to the membrane segment abutting cell i + 1: The number of auxin molecules in the cell changes as the result of their production, turnover and transport to and from the neighbouring cells. We assume that auxin is produced at a constant rate throughout the cell volume, which yields the propensity of event ∅ → A equal to Auxin turnover A → ∅ is treated as a random decay of auxin molecules, which yields propensity Consistent with Fick's law, the propensity of exporting auxin diffusively to a neighbouring cell is proportional to auxin concentration [a] and interface area s: A related formula describes the propensity of auxin transport facilitated by PIN, except that, in this case, the transport rate is modulated by the area density of PIN in the respective membrane: Analogously, propensity of transport to the right neighbour is Program 4a. Gillespie L-system productions implementing a stochastic model of auxin-driven leaf margin development.
Program 4b. Control of the production application in Program 4a.
The Petri net summarizing these processes is shown in Fig. 8, and the essential part of the resulting Gillespie L-system is given in Program 4a. The first rule in Gillespie group 1 captures the exocytosis and endocytosis of PINs, and the production and turnover of auxin. The application of this rule begins with a calculation of the number of PIN molecules in the cytoplasm using the equation where [PIN 0 ] is the volumetric concentration of PIN in the cell, assumed to be constant. The second Gillespie rule stochastically selects a transport event. This rule is centred on a wall rather than a cell, and requires slightly more complicated indexing than that in Equations (19-21 within that cell. This rule works in concert with the productions in group 2 to effect molecule transport between cells, in a manner similar to the diffusion-decay and reactiondiffusion Gillespie L-systems. The final element of the model is the development of the leaf margin, giving the simulated leaf its shape. Beginning with an initial shape resembling a leaf primordium, the cells elongate at a constant rate. Moreover, they are displaced in the normal direction at a rate proportional to the concentration of auxin, creating leaf lobes as described by Bilsborough et al. (2011) and Prusinkiewicz and Lane (2013). The rules implementing this growth are not shown in Program 4a, as they are not specific to Gillespie L-systems. However, cell division is described by the rule in group 3. Upon reaching threshold length x max the cell divides symmetrically. Auxin molecules are apportioned according to the size of child cells: due to the symmetry, each daughter cell inherits 1 2 of all molecules (only slightly more complex, but more consistent with the spirit of stochastic simulation would be to divide molecules A between left and right cell using a binomial distribution). Auxin concentration in the daughter cells is thus the same as it was in the mother cell. Likewise, Equation (22) apportions PIN molecules to the daughter cells proportionally to their volume. The numbers of PIN molecules allocated to the left and right membrane regions of each daughter cell are assumed to be the same as it was in the mother cell, so that the cell division does not disrupt the auxin flow. This implies preserving PINs in the existing membrane regions, and allocating PINs to the emerging wall between the daughter cells.
Growth and cell divisions could be simulated using stochastic productions as well, but the model achieves better performance by employing standard, deterministic L-system productions for this purpose. The resulting combination of stochastic and deterministic productions requires a careful scheduling of events (Lu et al. 2004). This is achieved by preceding the productions in Program 4a with a code controlling their execution, listed in Program 4b. Molecular processes are simulated first, by alternating between productions in ggroup 1 and group 2 as in Programs 2a and 3. This phase continues until the current simulation time t, returned by the GillespieTime () function, reaches time gt of the next growth-and-celldivision event. Productions in group 3 are then executed, and time gt is incremented by the predefined interval ∆gt. At this point, the simulation of molecular processes captured by productions in ggroup 1 usually resumes. It is possible, however, that simulation time t is still greater than growth-and-division time gt (this is particularly likely when the number of molecules in the model is small, implying large inter-reaction times and steps in the values of t). In this case, time gt is incremented again, and growth and cell division are simulated, until the condition t < gt becomes true. section S = 1 and were dividing upon reaching the threshold length x max = 22.5. Their initial lengths x i,0 were close to x max . The boundary cells were set to maintain zero auxin concentration, and the remaining six cells of the leaf primordium were initialized with 100x i,0 s auxin molecules each. Initially no PIN molecules were allocated to the cell membranes. Although these values-and details of the underlying equations-have been chosen arbitrarily, the model does illustrate the general applicability of Gillespie L-systems to the simulation of auxin-driven patterning processes, and their ability to capture random variations in this context.

Conclusions
We have proposed an integration of Gillespie's Stochastic Simulation Algorithm and L-systems as a method for simulating stochastic processes in structures with a constant or variable number of modules representing cells or higher-level compartments. While 'ordinary' L-system productions are applied to all modules in parallel, Gillespie-style productions are selected according to a set of propensity functions and applied to one module per simulation step. We have illustrated the operation of Gillespie L-systems with examples progressing from a single-compartment Lotka-Volterra model to diffusion-decay, reaction-diffusion and auxindriven morphogenetic processes. For simplicity, we have only considered linear structures (files of cells), although the formalism inherits from L-systems the capability of simulating branching structures as well. The combination of Gillespie's algorithm and L-systems makes it possible to account for the noise occurring in systems in which the number of molecules is small, and captures the variation in patterns and forms stemming from this noise. Prospective improvements and extensions include acceleration of simulations. One possibility is to limit the explicit computation of propensities to those affected by the previous simulation step, while propagating the remaining propensities intact (Gibson and Bruck 2000). The challenge is to automatically construct the dependency graph that would identify the propensities in need of updating, given an arbitrary Gillespie L-system. Other paths to acceleration are offered by improvements to the subvolume method (Elf and Ehrenberg 2004) and fast approximations of the Stochastic Simulation Algorithm (Gillespie 2007;Marquez-Lago and Burrage 2007;Lampoudi et al. 2009). Of interest is also an extension of Gillespie L-systems to two-and three-dimensional cell complexes (Desbrun et al. 2008;Lane 2015), which would allow for a stochastic simulation of processes taking place in growing tissues.

Sources of Funding
The support of this work by Discovery Grants 2014-05325 and 2019-06279 from the Natural Sciences and Engineering Research Council of Canada (P.P.), and the Plant Phenotyping and Imaging Research Centre/ Canada First Research Excellence Fund (P.P. and M.C.) is gratefully acknowledged.