VirtualLeaf: An Open-Source Framework for Cell-Based Modeling of Plant Tissue Growth and Development

• © 2011 American Society of Plant Biologists

Plant scientists have gathered a wealth of knowledge about plant growth and plant genetics. Yet, it is still impossible to predict by which mechanisms genetic changes affect plant growth and development. To unravel how genetic information determines the morphogenesis of plants, analyzing and reconstructing the dynamics of the genetic regulatory networks (Gonzalez et al., 2009) is only the first step (Merks and Glazier, 2005). To mechanistically predict morphogenesis from gene function, it is essential to have a detailed understanding not only of genetics but also of (1) how genetic networks regulate cell behaviors, (2) how cell behaviors lead to tissue growth and patterning, and (3) how the tissue-level phenomena, including spatial patterns (Swarup et al., 2005; Jönsson et al., 2006; Bayer et al., 2009) and strains and stresses (Green, 1999; Hamant et al., 2008), induce responses at the molecular level. How these multilevel feedbacks between levels of organization produce biological function and form is perhaps the most central question in systems biology (Noble, 2006).

Traditionally, biologists are familiar with presenting and developing hypotheses using “block-and-arrow” diagrams, symbolically representing identified and hypothetical gene interactions. However, despite their frequent use, it is not possible to understand the behavior of gene networks just by reading such diagrams, except for very straightforward cases. Mathematical models and computational models can bring the gene networks represented by block-and-arrow diagrams to “life” and can explain potential dynamic behaviors of the network even where precise quantitative information is unavailable. In this way, the biologist may discover that the hypothetical biochemical network cannot exhibit the expected behavior or that it occurs only for biologically unfeasible values of the reaction rates or binding constants. The researcher can then propose new interactions or system components and add them to the model to check if the revised network better reproduces the observations the original diagram was attempting to explain. Subsequent experiments will test the validity of the new elements in the hypothesis and can reveal a set of new observations that will be used to refine the mathematical model. This iterative approach of alternating experimental steps and model refinements is known as the “systems biology cycle” (Kitano, 2002). Mathematical and computational modeling supports experimental biologists in testing alternative hypotheses for “internal consistency” and helps them fine-tune their ideas.

A typical systems biology approach is taken in analyzing the de novo patterning of root hairs and in analyzing the regular arrangement of trichomes in the leaf epidermis of Arabidopsis (Arabidopsis thaliana). These two processes are controlled by the same molecular regulators (Schnittger et al., 1999; Pesch and Hülskamp, 2004; Kang et al., 2009). Based on a series of knockout experiments, an initial block-and-arrow model for the interaction of the transcription factors TTG, GL1, and TRY was proposed (Schnittger et al., 1999). By applying a local activation and lateral inhibition model for pattern formation, it was predicted that TRY is activated by GL1 and/or TTG and inhibits GL1 and/or TTG, two predictions that subsequently could be confirmed experimentally. Newly identified components have been added to the initial model (Pesch and Hülskamp, 2004), which was fine-tuned in a series of computer simulation studies (Benítez et al., 2007; Bouyer et al., 2008; Dupuy et al., 2008) that were validated experimentally (Kang et al., 2009). Ongoing cycles of computational modeling and experimental validation continue to fine-tune an awareness of epidermal patterning in Arabidopsis (Bouyer et al., 2008).

The analysis of epidermal patterning is one of several specific questions in plant morphogenesis that could be addressed with this approach, including gravitropism (Swarup et al., 2005), phyllotaxis (de Reuille et al., 2006; Jönsson et al., 2006; Smith et al., 2006), root tip specification (Grieneisen et al., 2007), leaf venation (Merks et al., 2007), mechanical signaling at the shoot apex (Hamant et al., 2008), lateral root initiation (Laskowski et al., 2008; Dun et al., 2009; Prusinkiewicz et al., 2009), and patterning of the cell wall (Tindemans et al., 2010). These projects demonstrate how plant biologists benefit from computational modeling by developing new ideas, insights, and hypotheses and how it helps predict missing system components. To date, modeling typically requires close collaboration between wet-lab scientists and computational biologists or mathematicians, even for testing relatively straightforward ideas or variants of existing models for which no new modeling and simulation codes and techniques need to be developed. Developments in bioinformatics demonstrate that, presumably, this only represents a transient state in the development of modeling. Initially, only computer scientists or computational biologists developed and applied sequence-based algorithms. Since then, experimental biologists have integrated several now-standard bioinformatics tools in their research, while computer scientists and mathematicians continue to develop new algorithms and methods. For example, to perform a BLAST search, experimental biologists no longer require collaborations with an expert bioinformatician. We hope that experimental biologists will soon also have a standard set of tools available for testing simple dynamic hypotheses on morphogenesis (Merks et al., 2006b). Apart from developing their own dynamic models and testing those in the laboratory, it will also help experimental biologists develop an intuition of how computational models can help them advance their research and assist them in framing a set of specific questions to be answered using more advanced models in collaboration with computational biologists.

A number of simulation systems exist for modeling multicellular plant systems. A popular and mature system is L-studio (Karwowski and Prusinkiewicz, 2004). As one of its features, L-studio implements vv-systems, a two-dimensional rewriting grammar to model cell division (Smith, 2006) derived from L-systems (Lindenmayer, 1968a, 1968b, 1975). L-studio and vv-systems have been applied in a number of recent studies on plant development (Smith et al., 2006; Bayer et al., 2009; Prusinkiewicz et al., 2009). L-studio has an intuitive user interface, and users can specify models using L-system or vv-system rules. In vv-systems, a morphological transformation of the tissue as a whole is specified by the user or by a model variable (Smith et al., 2006). The cell division algorithm then partitions the resulting space. This is a powerful approach if one is primarily interested in how growth can drive pattern formation (i.e. by increasing the domain in which pattern formation takes place; Crampin et al., 2002), such as leaf development (Runions et al., 2005), phyllotaxis (Jönsson et al., 2006; Smith et al., 2006), or trichome initiation (Dupuy et al., 2008) or if one is interested in how stable patterns can drive morphogenesis entirely downstream (see the emerging primordia in the phyllotaxis model of Smith et al., 2006). In contrast to such “morphostatic” development (Salazar-Ciudad et al., 2003), in many cases morphogenesis is a closed feedback loop of pattern formation, growth, and morphogenesis (Holloway, 2010); here, cell divisions guided by spatiotemporal patterns lead to domain growth that, in turn, changes pattern formation. In such “morphodynamic” systems (Salazar-Ciudad et al., 2003), the morphology emerges from a feedback loop between pattern formation and simultaneous morphogenesis, while chemical patterns drive morphogenesis by locally affecting cell behavior.

Other popular multicellular modeling techniques include the cell-based models (Merks and Glazier, 2005). Cell-based models derive morphogenesis from the genetically regulated behavior of individual cells that can be modeled as individual particles or as collections of particles. Chaste (for “cancer, heart, and soft-tissue environment”; Pitt-Francis et al., 2009) provides a set of libraries for developing biological simulations of animal tissues, including a cell-based simulation system for modeling tissue growth. It represents cells by its centers and connects cells with virtual springs. A disadvantage of such single-particle methods is that they cannot represent any morphological, subcellular detail, including cell polarity or cell shape. An example of a multiparticle method is the cellular Potts model (CPM; Graner and Glazier, 1992; Glazier and Graner, 1993). The CPM was developed originally for simulating animal development (for review, see Merks and Glazier, 2005) and was recently applied to plants (Grieneisen et al., 2007). In their model of root development, Grieneisen and coworkers (2007) propose that the auxin maximum at the root tip determines the location of the quiescent center in the root apical meristem. The auxin maximum results from a dynamic equilibrium between central auxin transport toward the root tip (acropetal transport), lateral auxin transport toward the shoot (basipetal transport), and auxin reflux from the lateral epidermis toward the stele at the root basal meristem. The same polar auxin transport mechanism also forms an auxin gradient along the root meristem, which is thought to regulate cell division, cell elongation, and cell maturation in the growing root via the PLT1 to PLT4 genes (Galinha et al., 2007). Although most of the work focuses on polar auxin transport on static cell lattices or cellular configurations derived from microscopic images, the CPM shows that the mechanism can produce realistic root division patterns while maintaining a dynamically stable auxin maximum at the root apex.

Because Chaste and the CPM in their present forms lack an accurate description of the physical properties of the cell wall, adjacent cells slide along one another, making their use in the modeling of growing plant tissues relatively limited. Therefore, most developmental phenomena simulated with Chaste and the CPM are examples of plastic morphogenesis, as it occurs in most animal tissues where the cells can migrate and move relative to one another. Because plant cells cannot migrate and do not slide along each other (an exception is the intercalary growth of phloem fibers [Ageeva et al., 2005]), plant morphogenesis depends exclusively on patterned cell division and cell expansion. This mode of development is often called symplastic growth (Priestley, 1930; Erickson, 1986).

Apart from “gluing” plant cells to each other and maintaining the shape of plant cells, cell walls determine the mechanical properties of plant tissues. The directed orientation of cell wall microfibrils determines the differential directionality of tensile strength in cell walls, allowing for anisotropic cell expansion that underlies organ morphogenesis (Green, 1980; Baskin et al., 2004). Moreover, cells respond to forces in the tissue, such as by setting their division plane parallel to the strain direction in the shoot apical meristem (Hamant et al., 2008). Thus, both the mechanical properties of the cell wall and the tight intercellular connections are crucial for understanding the morphogenesis of plant tissues.

A number of cell-based modeling methodologies consider detailed cell wall mechanics. CellModeller (http://www.archiroot.org.uk/doku.php/navigation/cellmodeller; Rudge and Haseloff, 2005; Dupuy et al., 2008) is a two-dimensional environment for plant tissue simulation. It has explicit representations of the cells, the cell walls, and the relative interactions between the cells and includes differential equation models of the dynamics of cellular properties, including biochemical networks. CellModeller describes cell walls as viscoelastic rods, which are maintained tensed by the turgor pressure. Cell wall dynamics is simulated with differential equations, while architectural, geometric models describe cell division and cell death. The model of Corson et al. (2009) also represents cell walls as viscoelastic beams while using an energy minimization approach to simulate cell wall dynamics and cell growth.

A shortcoming in most of the current modeling approaches is that they require a high level of mathematical and computer programming skills, which often keeps biologists from using them. Therefore, we introduce VirtualLeaf, a modeling environment that allows biologists to turn their hypotheses of multicellular development and patterning into a dynamic simulation model with little programming effort. For simplicity and computational efficiency, the underlying computational method describes cells in relatively little detail. Nevertheless, the model explicitly considers the balance between turgor pressure and cell wall tension. Our method relates to CellModeller and the method of Corson et al. (2009) in its representation of the cell walls and its consideration of turgor pressure. We illustrate the use of our simulation environment with a number of example simulations, including simple growing tissues and systems in which chemical pattern formation interacts with local growth.

MODEL DESCRIPTION

Cellular Mechanics

Our computational method represents cell walls in a tissue in two dimensions using polygonal finite elements (Fig. 1), an approach similar to that taken in other plant tissue models (Nagai and Honda, 2001; Rudge and Haseloff, 2005; Dupuy et al., 2008; Corson et al., 2009). A cell is bound by a set of elastic cell wall elements. To restrict relative cell movement, adjacent cells share a cell wall, where each cell wall consist of multiple cell wall elements, so that cell walls are flexible. A Monte Carlo-based energy minimization algorithm straightforwardly describes cell behaviors, including expansion, division, and active shape change. This approach conceptually derives from the CPM (Graner and Glazier, 1992; Glazier and Graner, 1993), which is most suitable for plastic, animal development (for review, see Merks and Glazier, 2005). We allow cell walls to move according to rules derived from the classic equations for turgor-driven cell expansion (Lockhart, 1965): (1) the intracellular turgor pressure exerts a uniform force on cell walls, attempting to enlarge cells; (2) the elastic cell walls counteract the turgor pressures; and (3) walls expand irreversibly if stretched over a threshold. Wall relaxation occurs under the influence of enzymes secreted into the wall that break specific bonds that cross-link cellulose microfibrils (Cosgrove, 2000) and is also called cell wall yielding (Lockhart, 1965; Cosgrove, 1985).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Cell-based modeling of plant tissue growth. Polygons represent cells, surrounded by cell walls, which consist of viscoelastic wall elements linked by nodes. A wall furthermore represents the membranes and transporter proteins on either side of the cell wall and the apoplast (i.e. the intercellular space). The algorithm iteratively displaces the nodes to balance the turgor pressure with the wall tension and possible extra constraints (e.g. a cell shape constraint). A, Overview of the simulation algorithm. The algorithm attempts to displace each node in random order; one such cycle is a Monte Carlo Step (MCS). MCSs are repeated until the energy fails to drop significantly. This iteration of MCSs is called a relaxation cycle. Once the system has relaxed, we simulate the biochemical networks and other cell behavior rules. These potentially affect cell properties (e.g. turgor pressure, wall extensibility, etc.), and a new relaxation cycle is required. B, Visualization of the simulation algorithm. Step 1, This move of a cell wall node relaxes local turgor pressure and minimally strains the cell walls; the move will be kept. Step 2, This move compresses the right cell and extends the left cell. The top wall element is overstrained, and the bottom wall element is compressed. Because this potential move would be energetically unfavorable, it will typically be rejected. Step 3, Each cell potentially contains a biochemical network described by a set of differential equations; the chemical species often regulate cell and wall attributes, including turgor pressure and wall extensibility. Step 4, A molecular species (shown in purple) channels back and forth between the cytoplasm and the cell membrane, as described by sets of differential equations. Step 5, Transport of the blue molecular species from cell to cell; diffusion passively transports molecules from high to low concentrations, and molecules at the cell membrane (e.g. PINs) actively transport molecules upstream or downstream.

Based on this set of rules, our algorithm defines for each cell i a resting area A_T(i) at which the cell’s turgor pressure balances the ambient pressure and a resting length L_T(j) for each wall element j; this is the length a section of cell wall would assume in the absence of turgor pressure. The actual cell area is a(i), and actual wall-element lengths are denoted l(j). We can describe the balance between turgor pressure and cell wall resistance in terms of a generalized potential energy or Hamiltonian (H):where indices i and j sum over all cells and polygon edges, respectively, λ_A is a parameter setting the cells’ resistance to compression or expansion, and λ_M is a spring constant.

The simulation uses Metropolis dynamics (Metropolis et al., 1953) to minimize the potential energy H: The algorithm iteratively selects a random node and attempts to move it in a random direction , where is a random vector chosen uniformly within the unit circle (i.e. ρ ∈[0, 1] and θ ∈[0, 2π]), and ξ is the step size that is fixed for a given simulation. The algorithm calculates the energy difference ΔH associated with the attempt and always accepts the move if it leads to an energy drop. To prevent the system from getting stuck in local energy minima, we also accept moves producing a small energy increase according to the Boltzmann probability function, P(ΔH) = e^−(ΔH/T). T is a parameter setting the amount of “noise” added in this way. To model wall relaxation (cell wall yielding), we introduce new nodes whenever a polygon edge’s length is stretched to the extent that it exceeds a threshold value ϑ_yield, typically 4 times the original length of the wall element. After the algorithm has attempted to move each of the nodes in the simulation in random order (a Monte Carlo step), it checks if the forces described in the Hamiltonian have balanced out sufficiently. Assuming that shape relaxation is fast compared with biological changes, the algorithm repeats this procedure until turgor pressures, cell expansion, and the resulting cell displacements have equilibrated sufficiently, and the energy difference drops below a relaxation threshold (ΔH < θ_H). After this relaxation cycle has completed, the algorithm solves the biologically motivated processes, including turgor pressure increases, cell divisions, active cell shape changes, and transport of chemical signals.

Thanks to the energy minimization formalism used for the model, we can include additional cellular properties in a straightforward way by describing them as energy constraints (e.g. body forces like gravity). For example, the length of cells (Merks et al., 2006a) can be constrained using:where H′ is the extended Hamiltonian, e(j) is the actual length of cell j, E(j) is the resting length of cell j, and λ_E gives the strength of the length constraint and can be interpreted as a deformability parameter. Assuming that cell shapes are ellipsoid, we estimate the cell lengths from the cells’ inertia tensors as with λ_b(j), the largest eigenvalue of cell j’s inertia tensor:where ∫∫D is the area integral over the cell and {x_c, y_c} is the cell's center of mass. Using Green’s theorem, the inertia tensor derives from the vertices of the polygon describing the cell. The algorithm reconstructs the inertia tensor from the raw moments of the cells’ coordinates, quantities that can be efficiently updated using local information after each move of the cell boundary. Thus, this algorithm produces a computationally efficient, instantaneous estimate of the cell lengths.

For cell division, we define a division axis and build a new cell wall (consisting of several wall elements) shared by the two daughter cells. The nodes of the parent cell are assigned to the two daughter cells. Typically, cell division is over the shortest axis, which is derived from the eigenvectors of the cell’s inertia tensor (Eq. 3).

Transport of Chemical Signals and Metabolites

Sets of ordinary differential equations in each cell describe the biochemical networks and genetic regulatory networks. Additional differential equations model diffusive transport across the two cell membranes and the cell wall separating adjacent cells. To model active transport by transport proteins (e.g. PINs), we define an extra object “wall” (Fig. 1B), which typically spans several “wall elements.” It stores the amounts of transport proteins present in the cell membranes at either side of the walls. Thus, the export rate of chemical A from cell 1 to an adjacent cell 2 would depend on the concentration of A in cell 1 and the amount of A’s export protein in cell 1’s membrane adjacent to cell 2. Although diffusion between apoplast compartments is not implemented yet, the wall elements can also describe an extracellular space (apoplast). For simplicity, we do not consider apoplasts in the examples described here. After each relaxation cycle (Fig. 1A), we dynamically construct a set of differential equations describing the intracellular reactions and transport of chemicals (with one equation for each chemical per cell and one equation for each chemical per cell wall) and numerically integrate them over a period Δt. This setup implies that the set of differential equations increases in size after each cell division, because division adds new cells and new cell walls. As a numerical integrator, we use a fifth-order, adaptive-step-size Runge-Kutta algorithm (Press et al., 1992).

Modeling Framework

VirtualLeaf separates the conceptual biological model from its implementation in terms of the algorithms for cellular mechanics and chemical transport as well as the visualization and steering of the simulations. The model definition itself is contained in a model plugin, a software module that is developed independently from the main framework. It is defined by combining software building blocks representing objects familiar to the experimental plant biologist, like cells, cell walls, molecules, transporter proteins, or signaling molecules. New model plugins are developed by adapting a small section of C++ code (Supplemental Text S1). The model framework runs fast (minutes rather than hours, depending on model complexity), enabling quick interactive desktop model evaluation. The models can also be run without a user interface, enabling high-throughput parameter screening on a computer cluster.

RESULTS

Based on the developed computational VirtualLeaf framework, it is possible to investigate various aspects of tissue growth. We will illustrate this with a number of example models that investigate emergent behavior of biological systems based on the characteristics attributed to individual component cells.

Example 1: Growth, Margin Stiffness, and Cell Shape

As a first study, we test the behavior of a simple model of plant tissue growth (Fig. 2A; Supplemental Video S1). The simulation was initiated with a single cell, whose target area A [and hence its “turgor pressure” (A_i − a_i)²] gradually increases at a constant rate. Cells divide over the short axis once their area a has doubled (for implementation details, see Supplemental Text S1, Tutorial 1). This basic simulation rule creates a rounded “callus” of plant tissue. In these simulations, we assumed that cell walls at the leaf perimeter were thicker and stiffer than the internal walls [λ_m(perimeter) = ζλ_m(internal)], with ζ = 2, the stiffness ratio between internal and perimeter walls. Perimeter walls are marked with blue nodes in Figures 1 and 2 and Supplemental Videos S1 to S3. This assumption agrees with the observation that margin walls in plant leaves are thick relative to walls of internal cells (Reinhardt et al., 2007) and that epidermal tissues are stiff relative to internal tissues in meristems. With these assumptions, our model rules produced a rounded aggregate of plant tissue. Interestingly, after reducing the stiffness of perimeter cell walls [λ_M(perimeter) = λ_M(internal)], a more “tumorous” morphology develops (Fig. 2B; Supplemental Video S2). In the absence of a neighboring cell or a thick cell wall resisting turgor, the perimeter cells “bulge” outward and expand faster than the internal tissue. Indeed, in plant leaves, the margin cells play a crucial role in morphogenesis (Reinhardt et al., 2007), which was proposed to be due to a biophysical function of the margin cells. In this simple model of plant tissue growth, we have not included the elongated shape of margin cells and we neglect many other mechanical features of plant cells (e.g. microtubules) and of the whole organ (e.g. the venature). Nevertheless, already in this simple model, the mechanical properties of the margin can regulate the shape of the tissue as a whole by constraining the expansion of the interior tissue. Figure 2C and Supplemental Video S3 illustrate the effect of a cell shape constraint (Eq. 2). In this example, the cells elongate into a direction that depends both on the initial shape of the cells and on the forces that the adjacent cells exert on them. Illustrating the shape constraint, this simulation is not intended to correspond to a real plant organ.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Symplastic cell-based plant tissue modeling. Cells adhere tightly to each other and expand at fixed rate. Cells divide once their area a_i has doubled. Simulation states are shown after 50, 100, 150, 200, 250, and 300 relaxation cycles (from left to right). A, Walls at the leaf perimeter have double stiffness relative to internal cell walls [λ_m(perimeter) = 2 λ_m(internal)]. B, If perimeter walls have “single stiffness” [λ_m(perimeter) = λ_m(internal)], perimeter cells bulge out, generating a tumorous morphology. C, Growing tissue with elongating cells. [See online article for color version of this figure.]

Example 2: Signaling Molecules Regulating Growth and Differentiation

Plant morphogenesis results from regulated cell divisions and pattern formation via polar or diffusive transport of chemical signals. Using the VirtualLeaf environment, it is possible to investigate the interaction between a diffusive morphogen and growth (Fig. 3; Supplemental Video S4). One of the cells produces a morphogen (e.g. cytokinin or auxin; concentrations indicated in shades of gray), which transports diffusively from cell to cell. If the concentration c of the morphogen is higher than a threshold, c > θ₁, where θ₁ is a concentration threshold, cells expand and divide once the area has doubled. At lower concentrations, θ₁ > c > θ₂, where θ₂ < θ₁ is a second concentration threshold, cells only expand. At concentrations c < θ₂, cells neither expand nor divide (i.e. they enter a maturation state). If the source cell divides, the source is assigned to one of the daughter cells at random. Note that in this illustration, we do not intend to model a biological system in particular; but even in this simple model, the cells naturally progress from a division phase to an expansion phase and finally enter a maturation phase, analogous to what happens in a root tip, for example (Beemster and Baskin, 1998). The model does not consider gravitropy; hence, the growth tip does not correct for the accidental sharp bend in Figure 3D. A local growth stagnation causes the bend: the morphogen concentration in the large expanding remains too high for the cell to enter division. This example simulation could form the basis for more complicated, realistic models focusing on the mechanisms responsible for morphogen gradient formation (Grieneisen et al., 2007) and on modeling the relation between cell division rate, cell expansion rate, and the size of the meristem and expansion zone (Chavarría-Krauser and Schurr, 2004; Chavarría-Krauser et al., 2005). Supplemental Text S1, Tutorial 2, shows how to construct such plant growth models involving signaling molecules with VirtualLeaf.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Diffusive morphogen stimulates growth. One tip cell produces a morphogen. At concentrations c > θ₁, cells expand and divide once the area has doubled, for θ₁ > c > θ₂, cells only expand, while for c < θ₂, cells neither expand nor divide (i.e. they go into a maturation state). Simulation state is shown after 100 (A), 500 (B), 1,000 (C), 1,500 (D), and 2,000 (E) relaxation cycles.[See online article for color version of this figure.]

A more realistic and complicated simulation involving chemical reactions that illustrates the use of the VirtualLeaf framework for developing new hypotheses of biological development involves studies of vascular patterning during Arabidopsis leaf development (Merks et al., 2007). The model predicts the existence of a traveling wave of auxin moving from the leaf tip to the leaf base. Auxin is concentrated into a peak by PIN1, which polarizes toward neighboring cells at a rate proportional to their concentration of auxin, a mechanism introduced by Jönsson et al. (2006) and Smith et al. (2006) for phyllotaxis. In our leaf venation model, auxin stimulates the production of auxin pumps. The pumps polarize toward the next cells and pump it onward. In this way, traveling waves form, which leave behind trails of polarized cells that we hypothesize to differentiate into the vascular system (Fig. 4; Supplemental Video S5). The long-term goal of this continuing project is to reconstruct in the VirtualLeaf environment how leaves develop in the interplay between patterned cell division and cell expansion, the transport of auxin and other signaling molecules between cells, and the physical characteristics of various cell types.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

A traveling wave of auxin forming a polar auxin channel. Auxin (shown in green) flows in from the leaf tip at the purple cell walls. PIN (red) localizes preferentially at cell walls adjacent to cells with the highest auxin concentration, as in recent phyllotaxis models. A, At 3 h, auxin accumulates at the leaf edge, and some auxin leaks into the lower cell layers. PIN1 production rises. B, At 13 h, auxin converges into a small peak, toward which the adjacent cells polarize. C, At 20 h, the auxin channel is completed, traced out by the traveling wave of auxin. See also Merks et al. (2007).

Example 3: Interplay between Growth, Cell Division, and Patterning

The next modeling example illustrates how tissue growth can interact with auxin-driven patterning. Here, each cell implements the upstream auxin transport model (Heisler and Jonsson, 2006; Jönsson et al., 2006; Sahlin et al., 2009) and lets local auxin concentrations drive cell expansion. For simplicity, we assume that auxin is produced along the tissue’s perimeter. Cells divide along their shortest axis when they reach twice their original size. After a small cell cluster has formed, auxin accumulates in the cluster’s center and induces further growth. New auxin maxima appear at the tissue perimeter once sufficient tissue has formed, inducing further, localized growth. This feedback between growth and pattern formation eventually forms a bulbous morphology (Fig. 5; Supplemental Video S6). Supplemental Text S1, Tutorials 3, 4, and 5, show in detail how such a model can be constructed in VirtualLeaf.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Interaction between auxin accumulation and growth. Auxin accumulates according to the model of Jönsson et al. (2006). Auxin induces cell expansion; cells divide over the shortest axis after doubling in size. Simulation state is shown after 10 (A), 30 (B), 50 (C), 70 (D), 90 (E), and 110 (F) simulation cycles.

In the previous sections, we illustrated the use of VirtualLeaf to model the biomechanics of plant tissue growth and its use in modeling pattern formation and morphogen-regulated growth. These relatively simple models give us more insight into interactions that exist between pattern formation, tissue growth, and tissue mechanics that drive the formation of plant growth. Supplemental Text S1 provides a step-by-step introduction to building models in VirtualLeaf.