Using cellzilla for plant growth simulations at the cellular level.

Cellzilla is a two-dimensional tissue simulation platform for plant modeling utilizing Cellerator arrows. Cellerator describes biochemical interactions with a simplified arrow-based notation; all interactions are input as reactions and are automatically translated to the appropriate differential equations using a computer algebra system. Cells are represented by a polygonal mesh of well-mixed compartments. Cell constituents can interact intercellularly via Cellerator reactions utilizing diffusion, transport, and action at a distance, as well as amongst themselves within a cell. The mesh data structure consists of vertices, edges (vertex pairs), and cells (and optional intercellular wall compartments) as ordered collections of edges. Simulations may be either static, in which cell constituents change with time but cell size and shape remain fixed; or dynamic, where cells can also grow. Growth is controlled by Hookean springs associated with each mesh edge and an outward pointing pressure force. Spring rest length grows at a rate proportional to the extension beyond equilibrium. Cell division occurs when a specified constituent (or cell mass) passes a (random, normally distributed) threshold. The orientation of new cell walls is determined either by Errera's rule, or by a potential model that weighs contributions due to equalizing daughter areas, minimizing wall length, alignment perpendicular to cell extension, and alignment perpendicular to actual growth direction.

Introduction

In recent years, there has been much interest in accurate multi-scale models of morphogenesis. Due to the various levels of complexity and the wide variety of tissue, there has necessarily been a trade-off between generality and specificity; an excellent review is given by Koumoutsakos et al. (2011). Of the most interest to plant biologists, perhaps, are platforms that describe the high-level structure of complete organisms, systems, or organs such as a meristem or sepal. Very few general purpose tools exist at all, and even fewer are specific to plants. L-System based tools, the best known being L-Studio (Karwowski and Prusinkiewicz, 2004), were among the first; they are ideally suited to branching structures because they are based on formal language theory (a grammar based on axioms, a short alphabet, strings derived from that alphabet based on specific production rules that are tuned to branching). L-system rules have been embedded in higher level languages such as C [specifically, cpfg, (Prusinkiewicz and Lindenmayer, 1990)] to allow the encoding of models and the description of geometrical and topological relationships. The Virtual Plant's OpenAlea Platform (Pradal et al., 2008) provides a general-purpose collection of simulation modules that use Python as the primary scripting language and allows L-systems to be incorporated (Boudin et al., 2012).

Outside of plant biology, much of the interest has focused on grid-based models such as the Cellular Potts Model (CPM) (Graner and Glazier, 1992), which discretize to the desired sub-cellular level of detail as collections of entities driven by particle-particle interactions; and finite element models (FEM) which instead discretize to a sufficiently fine web-work of straight lines under the influence of mechanical forces. Because of the molecular level of detail of these methods they are often primarily stochastic in nature (Gilllespie, 1976). There has been some effort to provide general purpose platforms; examples include MCell (Stiles and Bartol, 2001), designed originally to simulate the synaptic junction; Smoldyn (Andrews, 2012), primarily emphasizing nano-scale specificity; ChemCell (Plimpton and Slepoy, 2005), which models protein networks within cells; and CompuCell3D (Cickovski et al., 2007) which focuses on cellular function.

Of particular interest is the use of rule-based models, e.g., BioNetGen (Faeder et al., 2005, 2009) and NFSIM (Sneddon et al., 2011). In a rule-based model, molecules are considered structured objects and their interactions are described by rules for transforming these objects. Rule based models are interesting because every process that occurs within a biosystem, such as chemical and physical interactions and growth, development, and cell division and death, can be described by rules. Similarly, Dynamical Grammars (Mjolsness and Yosiphon, 2006; Yosiphon, 2009; Mjolsness, 2013) define models in terms of operator algebras of stochastic processes, with simulation algorithms derived from the composition and expansion of time-evolution operators.

At an intermediate level of complexity there are tools that describe tissue at the multicellular level, treating each cell as a well-mixed compartment. Cells can be described either as simple point (or spherical) objects connected by breakable springs (Jönsson et al., 2004, 2005a,b; Mjolsness, 2006), or at some level of geometric complexity with springs between polygonal vertices (Rudge and Haseloff, 2005; Sahlin and Johnsson, 2010), as is currently done by Cellzilla. The VirtualLeaf (Merks et al., 2011) provides an interesting hybrid that uses the CPM in combination with mechanical springs with a Markovian relaxation algorithm to describe tissue growth. Cellzilla has the distinction among these software implementations of extending the collection of Cellerator Arrows to include multicellular interactions (Shapiro et al., 2003).

Materials and methods

Cellerator reactions

Interactions are described in terms of Cellerator arrows (Shapiro et al., 2003); the basic canonical form is {X → Y, k}. When several species are interacting the arrow expression is written as where the e and f are stoichiometries. Each reactant is converted to a single term in a differential equation according to mass action kinetics: where e and f are the stoichiometries of species U in the reaction on the left and right hand side of the reaction, respectively, and k is the rate constant of the reaction. More complex expressions can be built from this canonical form to represent exact mass action descriptions of multi-species complex enzymatic reactions, as summarized in Table 1, and numerous regulatory approximations such as Michaelis-Menten, Hill functions, and MWC equations, as summarized in Table 2. Furthermore, a large number of exact enzymatic expansions have been implemented using the basic canonical form via the KMech toolbox, including BiBi, BiTer, BiUni, MulS, OrderedBiBi, OrderedBiUni, PingPong, PingPongTerTerF, PingPongTerTerR, RandomBiBi, TerBi, TerTer, UniBi, and UniUni reactions (Yang et al., 2005).

Table 1

Examples of Cellerator mass action arrow form expansion, from the base from .

Arrow form*	Expansion
e1X1 + e2X2 + … ⇄	e1X1 + e2X2 + … → f1Y1 + f2Y2 + …
f1Y1 + f2Y2 + …	f1Y1 + f2Y2 + … → e1X1 + e2X2 + …
X⇄ℰY	X + ℰ ⇄ X_ℰ
	X_ℰ ⇄ Y + ℰ
X⇄ℰℱY	X⇄ℰY and Y⇄ℱX
X⇌ℰY	X + ℰ ⇄ X_ℰ
	X_ℰ ⇄ Y_ℰ
	Y_ℰ ⇄ Y + ℰ
X⇄Y⇄Z⇄⋯ℰ	X⇄ℰY, Y⇄ℰZ,…
X⇄Y⇄Z⋯ℰℱ	X⇄ℰℱY, Y⇄ℰℱZ,…

The complete syntax of an arrow form is {arrowform, k.

Table 2

Cellerator user defined, regulatory, and enzymatic arrow forms.

Arrow forma,b	Typical ODE Termc
{e · X ⇒ f · Y, g(X)}	[Xi]′ = −eig(X), [Yi]′ = fig(X)
{X↦ℰY, Hill[v, n, K, a, T]}	[Y]^′=v[ℰ] (a + T · [X])^nK^n + (a + T · [X])^n
{X↦ℰY, GRN[v, T, n, h]}	[Y]^′=v[ℰ]1+exp(−h−T · [X]^n)
{X↦ℰY, SSystem[τ, C_+, C_−, n^+, n^−]}	[Y]^′=[ℰ]τ{C_+∏_iXini+−C_−∏_iXini−}
{X↦ℰY, NHCA[v, {T^+, T^−}, n, m, k]}	[Y]^′=v[ℰ]∏_i(1+Ti+[X]in_i)^mk∏_i(1+Ti−[X]in_i)^m+∏_i(1+Ti+[X]in_i)^m
{{X, Y} ⇒ rational[a, d, m, n]}	[Y]^′d=a_0 + rest(a) · [X]^nd_0 + rest(d) · [X]^n
{X ↦ Y, USER[v, T, n, h, f]}	[Y]′ = vf(h − T · [X]n)
{X⇒ℰY, MM[K, v]}	[Y]^′ = −[X]^′ = v[ℰ][X]K + [X]
{X⇒ℰY, MWC[k, n, c, L, K]}	[Y]^′e=−[X]^′= k[ℰ]α(1+α)^n−1+Lαc(1+αc)^n −1(1+α)^n+L(1+αc)^n−1

.

.

[ℰ] with 1.

rational, rest(.

.

.

.

.

Cellzilla is implemented in Mathematica and is invoked using the standard notebook interface. The details of numerical integration are normally hidden from the user, e.g., models (such as (30)–(32)) are submitted to the kernel by the Cellzilla Grow command and Cellerator run along with a list of simulation control parameters and initial conditions. These front-end commands directly invoke Mathematica's NDSolve. Normally, the default solver in NDSolve is chosen, unless the user selects otherwise; any control option for NDSolve may be passed along by Grow command. For ordinary differential equations, NDSolve switches between a non-stiff Adams method and a stiff Gear Backward Differentiation formula. However, users may optionally change the parameters or choose another solver.

Tissue description

Tissues are described by a polygonal lattice, with each lattice cell representing one biological cell. The Tissue data structure consists of: (1) a vertex list V, where each vertex V is an (x, y) pair; (2) an edge list E, where each edge E = (i, j) is a pair of integers giving indices of vertices at the endpoints of the edge; and (3) a list of cells C, where each cell C = {k1, k2, k3, …} is an ordered list of edge indices. Externally the tissues may be saved (either read or written) as CSV files, either as lists of vertices, edges, and cells, or in a flattened versions with the edges omitted and the cells represented as ordered sequences of vertices. Internally the edges are always reconstructed since it is more efficient computationally to always have the edge information available, but it is user taste which I/O format is used. It was decided to use the CSV format for these files because of the wide availability of parsing tools and the ease of human readability should that be necessary.

Species in different cells are referenced by an index; e.g., the reaction X[17] → Y[17] takes place in cell 17. When a reaction network is expanded in every cell in the system it is not necessary to repeat this manually, as this is done automatically. Constituents in different cells can interact in the following ways: (a) diffusion; (b) action at a distance; (c) transport across the cell wall. Each of these may be specified in the model by one of the additional Cellzilla Arrow forms listed in Table 3. Such interactions are allowed to depend on a function f(i, j, k) that depends on the properties of the constituents of the cells (i = present cell number, j = connecting cell numbers) and cell wall (k) between cells i and j. The basic single cell models may be input and output as SBML files using the Cellerator/MathSBML extensions for SBML (Hucka et al., 2003; Shapiro et al., 2004).

Table 3

Additional Cellzilla arrow forms not recognized by Cellerator.

Arrow form*	Description
{X → X, Diffusion[PI, PO]}	Diffusion of X through the tissue. PI is the permeability of internal cell walls; the optional PO is the permeability of tissue boundary cell walls. Each may be specified as f[i, j, k] where i, j, k are cell and wall indices.
{X → X, Transport[fout, fin]}	Controlled transport of X across the cell wall.
{X ↦ Y, IGRN[v, T, n, h]}	[Y[j]]^′=v1+exp(−h−T · [X[i]]^n)
{cell → cell, Grow[…]}	Specification of cell growth parameters: Pressure, growth rate, spring constant.
{cell → cell + cell, model[…]}	Specification of cell division model and parameters.

Except for the IGRN the arrows here are longer versions of the right-pointing arrows used by the canonical Cellerator mass-action expansion.

Diffusion across cell boundaries is implemented according to Fick's law, so that the flux J through any membrane (e.g., in molecules/(cm2-s)) with diffusion constant D (e.g., in cm2/s) is

Defining the membrane permeability as β = D/δ (e.g., in cm/s), where δ the membrane (or wall) thickness gives where ΔX is the concentration difference across the membrane. Let cell i have area A and depth d (orthogonal to the simulation); the the volume of cell j is V = Ad, and the area of the cell wall between cell i and cell j is ℓd where ℓ is the length of the wall between the cells. The flux across wall k into cell i is (for constant area):

Therefore

Rather than implementing this equation directly, it is implemented as equivalent Cellerator reactions where f is the input concentration-dependent permeability. Diffusion is specified to Cellzilla by incorporating arrows of the following form into the model: where f is either a Mathematica pure function or a function that has been defined previously in the simulation.

Action at a distance is technique for describing the way constituents in cell i can affect constituents in another cell j without specifying any of the intermediate reactions. We restrict action at a distance to adjacent cells with Cellerator GRN-like reactions written as {X ↦ Y, IGRN[v, β, n, h]} which means that constituent X affects constituent Y in neighboring cell j according to (Mjolsness et al., 1991)

Facilitated membrane transport from cell j to cell i is described by the equations where f(i, j, k) is the positive outward molecular flux from i through edge k, to j, and f(i, j, k) is the positive inward molecular flux to cell i through edge k, originating from cell j. In most cases one will only want to specify one of f or f. These are implemented internally via the Cellerator reactions

Transport reactions of this sort are specified to Cellzilla by including arrows of the form where the function fout should be set to zero if only f is utilized.

Growth model

Cellzilla implements two types of time-dependent simulations: static, and growing. In static simulations the shape of the tissue and its component parts do not change but its constituents are allowed to vary as described previously. In a growing tissue, the shape of the cells are also allowed to evolve with time. Cell growth is described by associating a Hooke's law spring potential of the form with the edge connecting vertices x and x. Here δ is an equilibrium length assigned to the edge, ℓ is the actual length, and k is a constant. When the wall is under compression (so that δ > ℓ) there will be a force, acting along the length of the edge, pushing the two vertices apart; when the wall is extended (δ < ℓ), the force will tend to pull the vertices toward one-another. The magnitude of this force is equal to the negative gradient of V. In addition, a pressure P associated with each cell a is, is applied outward at each vertex. The net force on each vertex is proportional to each wall incident on that vertex and, and is split evenly between the vertices. Then equation of motion for vertex x is then where is a unit vector pointing from x to x. The sum in the first term is over all the neighbors j of vertex i. In the second term, n is an outward pointing unit-normal vector from cell a, normal to edge ℓ, and the sum is over all neighbors j of i and and over all cells a in incident at vertex i. The pressure force will cause the springs to extend, simulating cell growth. The resting length is allowed to increase linearly at a rate proportional to the extension beyond resting length, where

If the spring is not extended, then growth does not occur. Equations (15) and (16) capture the phenomenological behavior observed in nature by both plant and animal cells. Pressure drives cell expansion; as pressure increases, cells expand more quickly. The growth rate is limited by the spring force. The dynamics of (15) can be solved exactly for a square cell to give where L is the cell perimeter and Δ=∑δ is the sum of the resting lengths. Consequently, when μ = 0 in (16), the sigmoidal growth pattern of plant cells is observed [e.g., when extensibility decreases linearly over time and osmotic pressure is held constant, as in Figure 1 of Lockhart (1965)]. Comparing (18) to the Lockhart equation ℓ′/ℓ = Φ(P − P) we see that the k and δ are parameters that can tuned to fit effective extensibility Φ and yield pressure P; for constant P and k sigmoidal behavior can be tuned for P < 2k in square cells. Additionally, allowing μ > 0, along with the spring force produces a more general growth model that is more generally applicable, not just in plant tissue, as the springs can be cut beyond a specified threshold, thereby removing cell-cell interactions (Shapiro and Mjolsness, 2001).

Let r be the index of the edge connecting vertices i and j, so that δ and ℓ are short notations δ and ℓ. The spring dynamics of (16) are implemented internally as where is the pth Cartesian component (p = 1, 2 for x or y) of vertex x, and k[r] is a spring constant whose value may depend on the properties (e.g., constituents of) the cells abutting edge r. Similarly, the growth of edge δ described by (16) is implemented internally as the Cellerator reactions where p and q are the indices of the cells that about edge r. The pressure force in (15) on vertex x due to cell a on the edge connecting vertices x and x is implemented internally as where k = 1, 2 (to indicate x or y Cartesian component); P[a] is the pressure in cell a; and n[a, i, j, k] is the k component of a unit normal vector to edge ℓ pointing outwards from cell a.

Cell growth is specified in a Cellzilla model by a reaction of the form where the arguments to Grow specify growth parameters such as the dependence of k, P, and μ on cell constituents (see Table 3). Chemical concentrations change during growth in each cell occur even in the absence of reactions. If there are n molecules of X in volume V then [X]′ = (n/V)′ = (Vn′ − nV′)/V2 = n′/V − [X]V′/V. The second term gives the correction in [X]' due to volume changes.

Cell division

Division occurs when a cell's area passes a threshold. Upon birth, each cell is assigned a threshold that is distributed normally (with the mean μ and standard deviation σ as optional control parameters). Chemical concentrations are distributed equally between the child cells (so that the chemical amounts are proportional to cell area). A single (linear) near cell wall is placed according to one of two user-selectable model: the standard modern interpretation of (Errera, 1888) and a potential model. In the modern interpretation of Errera's rule, the shortest wall that divides the two cells in half (by area) is chosen. [Technically, this is not Errera's rule, which only defines the shape of the cell wall, once the endpoints are already know; however, the area-equalization constraint is typically added to provide this boundary condition. (Smith, 2001; Besson and Dumais, 2011; Prusinkiewicz and Runions, 2012)] In the potential model (2010, 2010) a function is minimized over the central angles θ1 and θ2. These give the central angles of the end points of the new wall measured from the cell centroid. Here w is a weight vector, and V represents each contributor to cell division, where i ∈ {A, L, e, g}, as described in the below.

The area potential V is minimized when the cell divides in half; if the daughter cells have areas A1 and A2, respectively, then we define

The function is squared to improve computational stability near the minimum (which would otherwise have a non-differentiable corner there). The disadvantage of this potential is that it does not have a unique global minimum, i.e., any line of cell division will that divides the area in half will give a value of zero. Thus, the area potential must be tempered by either an additional potential function (such as the perpendicularity and/or length potential) or an additional heuristic to select the desired minimum value that does not require a unique minimum (e.g., randomly select division direction from amongst all equivalent minima).

The length potential V will be minimized when the new cell wall is most closely aligned with the shortest possible diameter d that this, the shortest line segment dividing the cell passing through the cell center. If d is the length of the new cell wall, then V is given by where Δ is the shortest distance between the new wall and the cell center, and ϵ is a tunable parameter.

For V and V we define a perpendicularity potential V⊥(v) that is minimized when the new wall is perpendicular to particular unit vector v. Let W be a unit vector parallel to the new wall. Then where ϵ⊥ is a parameter. Letting e be the direction of maximal cell extension and g be the direction of maximal cell growth, we then define so that V and V are minimized when W is most nearly perpendicular to the directions of maximal extension and maximal cell growth, respectively. The direction of maximal extension is taken as the unit eigenvector corresponding to the larger eigenvalue of the covariance matrix M = (XX)/(n − 1) where n is the number of cell vertices; and X = [x − x|y − y], where x and y are column vectors of cell vertex coordinates {(x1, y1), …, (x, y)} and x and y is their mean. The direction of instantaneous maximal growth is found in the same manner, using the velocities of the vertices. The covariance matrix M is calculated using the Mathematica function Covariance.

Cell division is specified in a Cellzilla model with the arrow where model is either ErreraModel or Potential, and the arguments specify the threshold variable, mean, standard deviation, and weight vector (for the potential model). A more generalized version of this notation has been introduced by Yosiphon (2009).

Results

Templates and the brusselator

Cellzilla has variety of shapes that can be used for basic simulations such as rectangular and hexagonal arrays, as well as circular, semicircular and parabolic templates that can be populated with randomly placed Voronoi centers. Alternatively the user can supply a template of his or her own consisting either of cell enters (in which case the walls will be interpolated with a Voronoi algorithm) or cell walls as described previously. Here we present the use of several of these templates (hexagonal, Voronoi, and user-supplied) to implement a common reaction-diffusion system.

The Brusselator (Prigogine and Lefever, 1968) is frequently cited in mathematical modeling because it consists of a system of chemical reactions that in the appropriate parameter regime will maintain sustained oscillations. When combined with diffusion such a system can also be used to establish a wide variety of interesting patterns such as stripes, spirals, and central maxima. The establishment of these different patterns depends on the choice of system geometry, boundary conditions, and parameter values. A diffusible Brusselator is easily implemented in Cellzilla with where a, β, b, c, D and D are tunable parameters. As illustrated in Figure 1, the ratio of the diffusion constants will change the number of maxima achieved. In the second row of the figure we see that the geometry is also significant. While the qualitative features of each of the results D-F are identical, their symmetry becomes more and more broken as the symmetry of the template becomes lost. With the exception of the diffusion constants, all other parameters were identical through each of these simulations. We are particularly interested in the parameter set shown in Figure 1A, because it can be used as described in the following section to establish an organizing center for simulation of the WUS/CLV network.

Figure 1

Results of Brusselator simulation showing affect of ratio of diffusion constant and geometry on simulation. (A–F) The concentration of species A in the brusselator (30) is shown, with higher concentration given in dark blue, and zero concentration in white. (A) D; (B) D, DB = 0.1; (C) D, D; (D–F) all use D = 0, D = 0.5 on different cellular teimplates. (G) Time course of the concentration of species A for all cells in simulation B (total of 199 cells). Each curve gives the concentration for different cell. Simulation E uses a Voronoi template of 500 randombly placed centers, and F uses a actual Arabodopsis meristem L1 segmentation. Parameters: a = 0.1, β = 0.1, c = 0.1.

Establishment of stem cell niche

In Jönsson et al. (2005a,b) we presented a predictive model of feedback interaction between the Wuschel (WUS) and Clavata3 (CLV) signals in the shoot apical meristem. This simplified model was able both to organize the WUS expression domain and to predict the reorganization due to the removal of the CLV signal from the WUS domain as seen in experiments when cells are ablated. This model uses a reaction-diffusion mechanism to induce WUS; the pattern is induced by a Brusselator. The original model relies on a diffusible parameter Y that is produced only in the L1 layer of a slice. We present an implementation in which our slice has the L1 layer omitted, and replace this with a boundary condition in which Y is held fixed, and allowed to diffuse inward. Assuming the Brusselator is implemented by (30), the Cellzilla network for the WUS activator is given by where v, T, T, h, k, k and D are tunable parameters, and the control word sigma tells Cellerator to replace the usual logistic control function f(x) = 1/(1 + e−) with (in fact, any monotonic increasing saturating function would work). The results illustrated in Figure 2 show that both the original central maximum (in wild type) and dual, smaller maxima result (in the ablation experiment) as modeled previously. In addition, we show the steady state distribution of the constituent Y in Figure 2C, illustrating how it forms a ring in the outer cells abutting L1 and decreasing inward, as desired.

Figure 2

Results of Wuschel simulation using Cellzilla. Species concentration is illustrated at the end of the simulation time course. The concentration of Wuschel is shown in blue, and protein Y is shown in green. Darker colors indicate higher concentrations, and white indicates a zero concentration. Steady state concentration of (A) Wuschel in wild-type simulation; (B) Wuschel in ablated meristem; (C) signal protein Y in ablated meristem. Parameters: v = 0.1, T = −25, T = 0.5, h = 0, k = 0.1, k = 0.1, d = 0.5, D = 2, D = 1.5, D = 15.

Growth induced by organizing center

The Brusselator is a useful mathematical/computational artifice for establishing patterns that can be otherwise studied but its biological meaning becomes lost if the variables in the equations do not have biological analogues that are present in the actual tissue. It is more meaningful if the organizing tissue can be established based on specific networks whose constituents have been observed and whose interactions are believed to be present, although this may at times be more computationally intensive. For example, Nikolaev (Nikolaev et al., 2007, 2013) has shown that a combination of reaction-diffusion and feedback in the WUS/CLV network is sufficient to establish a stem cell niche. In a one dimensional dynamic model including cell division, Chickarmane et al. (2012) has shown that negative feedback between WUS and cytokinin synthesis may be sufficient for maintenance of this niche as the tissue grows. Here we present for illustrative purposes of Cellzilla capability a simplified version of a Chickarmane-inspired model, in which two diffusible species are used to establish the pattern: U (e.g., that may be part of the cytokinin network), which is produced only in the cells at the tip of the meristem; and a second species V that is produced in the L1 layer (e.g., that may produced as part of the CLV network). The species that represents the organizing center is called W; U and V then activate and repress W, respectively, while W is self-activating, perhaps through an intermediate. Positive feedback of W onto V is provided by a third diffusible species X; and the epidermis is impermeable to U, V and W. Finally, the constitutive degradation of W is slightly enhanced in the L2 layer, but occurs everywhere. The network is where k, k, k, k, k, k, k, k, k, h, h, v, v, T, T, T, T, T, D, D, and D are tunable parameters; L1[t], L2[t], and Tip[t] are built in indicator functions for these cell locations; and the functions f and fμ describe pressure and growth feedback; e.g., P[i] = p0 + p1W[i] and μ[j] = μ0 + μ1(W[i] + W[j]) (where p0, p1, and μ0 are tunable control parameters). The user would type these function definitions either in the lines before the model or in line as lambda functions in place of the function references in the model. Simulation results are shown in Figure 3. As seen in Figures 3A–C, the stem cell niche is maintained through at least 500 cell divisions.

Figure 3

Maintenance of organizing center during growth simulation. (A) Initial distribution of W; (B) Distribution after 250 cell divisions; and (C) Distribution after 500 cell divisions. (D) Cell lineages; each color shows a different clonal population corresponding to all descendants of a particular cell in (A). Parameters: k = 2, k = 0.2, k = 1, k = 0.25, k = 1, k = 0.05, k = 1.5, k = 0.1, D = 10, D = 0.5, D = 0.5, v = 1, v = 1, v = 1, h = 0, h = 0, h = 0, T = 4, T = 4, T = 22.5, T = −25, T = 27.5, p = 0.001, p = 0.004, μ0 = 5 × 10−6, μ1 = .004.

Discussion

We have illustrated that meaningful quantitative results for plant morphodynamics can be obtained using a simple polygonal tissue model coupled with a spring growth equations. In particular, our implementation utilizes and extends an existing arrow-based computational framework that is easy and intuitive to use. This framework is built within a standardized computer algebra system (Mathematica) that is widely available and provides access to a wide selection of analytical tools. In addition we believe that our framework is generalizable and extensible to the wider world or rule-based systems. While more detailed particle-based or molecular dynamics frameworks will certainly produce more accurate frameworks their ultimate extensibility is limited by CPU availability and time constraints. Our implementation provides a useful platform for rapid model development and testing that can easily by transformed to one of the more detailed frameworks, if desired, once suitable results are obtained.

Cellzilla can be used for 2D simulations of plant tissues at the multicellular level. Because models can be rapidly built and tested it allows one to quickly test developmental models both in steady-state and during growth. However, several improvements are planned for future versions. In the current version, only simple transport and diffusion are considered; however, in many cell types, not just plant cells, osmotic and electrical gradients are significant. We plan to implement rules to incorporate both features in future releases. Furthermore, with the addition of osmotic models and pressure gradients we plan to add additional plant-specific growth models (Lockhart, 1965; Ortega, 1985; Geitmann and Ortega, 2009). These can be optionally used in place of the phenomenologically-based spring-based growth model. Additions to the cell division model will be included. For example, the constraints used in the Errera implementation of linear walls and equal areas can be relaxed to quadratic and circular arcs, and the areas can be randomized in a cloud about equality.

A three-dimensional model is also planned, although we expect that this will have significantly greater computational demands. To avoid mathematically unstable solutions in three dimensions the spring model will need to be modified or replaced, most likely with an elastic dynamics model incorporating pressure and stress tensors. An alternative is the use of triangular springs (Delingette, 2008).

Data sharing

All of the software described here is open source Mathematica (GPL license) and freely downloadable from launchpad at https://launchpad.net/cellerator. The software is fully documented and all examples are available at the project website https://www.cellzilla.info.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.