Shaping the Organ: A Biologist Guide to Quantitative Models of Plant Morphogenesis

Organ morphogenesis is the process of shape acquisition initiated with a small reservoir of undifferentiated cells. In plants, morphogenesis is a complex endeavor that comprises a large number of interacting elements, including mechanical stimuli, biochemical signaling, and genetic prerequisites. Because of the large body of data being produced by modern laboratories, solving this complexity requires the application of computational techniques and analyses. In the last two decades, computational models combined with wet-lab experiments have advanced our understanding of plant organ morphogenesis. Here, we provide a comprehensive review of the most important achievements in the field of computational plant morphodynamics. We present a brief history from the earliest attempts to describe plant forms using algorithmic pattern generation to the evolution of quantitative cell-based models fueled by increasing computational power. We then provide an overview of the most common types of “digital plant” paradigms, and demonstrate how models benefit from diverse techniques used to describe cell growth mechanics. Finally, we highlight the development of computational frameworks designed to resolve organ shape complexity through integration of mechanical, biochemical, and genetic cues into a quantitative standardized and user-friendly environment.


INTRODUCTION
Plant development has inspired the interest of scientific minds since antiquity. The first attempt to formulate plant growth into a mathematically coherent framework was given by Thompson (1917) in his landmark book On Growth and Form. Idea of Thompson (1917) was that morphogenesis could be summarized as a series of coherent geometrical transformations leading to the spacious diversity of biological forms. The concept of morphogenesis is therefore quite general; we usually define morphogenesis as a recipe to build an organism with elements such as individual cells, genes products and biochemical signals. While many cells proliferate to recreate the organism's adult shape some of them may differentiate into specialized tissues. Understanding the rules behind this decision-making process is at the core of organ patterning. To grasp the principles of morphogenesis, one has to consider different scales of organization such as growth mechanics, biochemical reactions, and genetic blueprint. Given the sheer number of elements involved in a morphogenetic pattern formation (i.e., the plant embryo may contain tenths of cells), it is impractical to analyze it without the support from computational tools.
In biology, computer models of patterning made their appearance in the 1960s with the widespread use of computer algorithms. Most of these models were based on static templates and did not include cell growth dynamics (Turing, 1952). Remarkably, Ulam (1962) showed that incredibly complicated forms and structures could be generated using cellular automata. Over the years, many researchers have recognized the importance of positional information in morphogenesis, as growth and cell division affect chemical gradients by diluting and degrading biological molecules (Wolpert, 1969). One of the earliest computational models of expanding plant tissue was proposed by Korn (1969). This model simulates the expansion of the alga Coleochaete scutata by applying a series of rules for cell growth and division over a hexagonal lattice (Korn, 1969). A similar work on two-dimensional growth of cell populations was previously introduced by Eden (1961). Pattern generation was a popular topic in computer science in that period and researchers were eager to find applications in biology and other fields. Cohen (1967) published a computer program capable of producing branching patterns reminiscent of leaf vascularization. These initial attempts to represent descriptive growth were not formalized into a comprehensive computational framework until the introduction of L-systems by Lindenmayer (1968). L-systems quickly gained popularity in the 1970s and until today are one of the favored methodologies for modeling plant architectures. Similar approaches to L-systems were developed for models of tree-like branching formation (Honda, 1971(Honda, , 1983. Simultaneously, early applications of continuous mechanics methods [i.e., finite element method (FEM)] were applied to green algae development (Niklas, 1977). These were followed by the implementation of anisotropic deformation and stress-strain relation in models of plant growth (Silk and Erickson, 1979;Hejnowicz and Romberger, 1984;Peters and Tomos, 1996). Furthermore, the relentless increase of computational power of modern machines allowed the definition of increasingly complex structures, such as tissues and entire organs (Glazier and Graner, 1993).
These initial successes in modeling plant development motivated the incorporation of biomechanics properties to further increase the realism of plant form generation. Lockhart (1965) proposed a model to predict cell wall growth rate from internal turgor pressure by formulating a set of rules regulating the behavior of an idealized cell wall, which will later become the de facto standard for mechanical modeling of plant cell elongation. However, these formulas could not represent important mechanical aspects exhibited by living cells such as water transpiration, wall stress relaxation, pressure relaxation, and elastic deformations (Geitmann and Ortega, 2009), which were addressed by follow-up studies: Cosgrove (1986) extended Lockhart's paradigm (which only applied to single cells in isolation) to multicellular organization; Silk and Wagner (1980) proposed a non-compartmented continuum model; Ortega (1985) augmented Lockhart's equation with an elastic component, while Veytsman and Cosgrove (1998) rederived the formula in terms of thermodynamics of polymer networks. These studies greatly improved Lockhart's initial idea and transformed it into a more flexible framework suited for complex environments.
In the last decade, we have witnessed an exponential growth of theoretical and experimental studies incorporating either biochemical (Leyser, 2018) or biomechanical ) aspects of plant organ growth. The continuous feedback between these two signals can potentially give rise to outcomes hardly predictable through simple human intuition, and thus safeguarding the necessity for even more advanced computational modeling techniques. To guarantee reproducibility and standardization, several pre-packaged and ready-to-use modeling frameworks have been developed (Prusinkiewicz et al., 2000;Pradal et al., 2008;Merks et al., 2011) which provide an interactive environment for biologists lacking programming knowledge.
In this review, we evaluate common biological issues and bottlenecks in modeling plant organ forms as well as their implications with respect to model realism. Finally, we present some of the most popular modeling frameworks that have been developed in the attempt to solve these issues.

GRASPING ON THE COMPLEXITY OF PLANT ORGAN STRUCTURE
In general, computational models are designed to serve two main purposes: describing a natural phenomenon of interest to gain insight into the mechanics of the process or to make predictive statements about an idealized hypothesis to forecast the result of in vivo experiments. There are several recurring questions in plant growth and morphogenesis that have been under the scrutiny of researchers for a long time and that could be now addressed with quantitative computer models.
Plants possess apical/basal polarity axes (Jürgens, 2001). Single cells can expand either isotropically or anisotropically ( Figure 1A), and cell polarity contributes to the collective choice of single cells which eventually determine the future shape of the plant (Wabnik et al., 2013). Some plant organs such as the root, root hairs, fruit, stems, and pollen tubes display a clearly anisotropic shape (Baskin, 2005). Growth anisotropy correlates with the direction of cellulose microfibrils deposited during cell wall synthesis (Li et al., 2014), but the exact mechanisms of how anisotropic growth is achieved are largely unknown. Cellulose fibers are very stiff and can be commonly found disordered or oriented along a preferred cell axis, as a result the cell wall remains stiff along a specific direction (Cosgrove, 2005). Microfibrils are often parallel to cortical microtubule orientation and it is thought that cellulose microfibril synthesis is guided by cortical microtubule tracks (Paradez et al., 2006). It has been shown that cortical microtubules align with maximal tensile stress in plant tissues (Hamant et al., 2019), indicating the existence of a feedback loop which reinforces the cell anisotropy against principal stress directions ( Figure 1A). Surprisingly, it has been suggested that even isotropic growth could generate anisotropic patterns as a result Frontiers in Plant Science | www.frontiersin.org 3 October 2021 | Volume 12 | Article 746183 of differential growth rates between adjacent tissues (Kennaway et al., 2011;Marconi et al., 2021). Other studies have also hypothesized the existence of morphogen-driven growth polarity fields (Mansfield et al., 2018). The factors partially involved in anisotropic growth are multiple and are not limited to cytoskeleton configuration. Cell elongation is mainly driven by two mechanisms: cell wall deformation by internal turgor pressure and cell wall growth mediated by enzymatic reactions [involving expansins, xyloglucan endotransglucosylase/hydrolase (XTH), and pectin-modifying enzymes (PMEs)] and favored by hormones accumulation (like auxin, gibberellins and abscisic acid; Cosgrove, 2016;Smithers et al., 2019). The latter is commonly known as the acid growth hypothesis ( Figure 1B; Rayle and Cleland, 1992). Turgor pressure is created by the continuous uptake of water from the external environment, which in turn exerts physical stresses on the cell wall (Cosgrove, 1993). The cell wall is incredibly rigid, capable of sustaining extremely high internal pressure (up to three atmospheres), and it is generally regarded as a viscoelastic material (Geitmann and Ortega, 2009). For the sake of simplicity, turgor pressure can be considered constant during cell growth (Schopfer, 2006), and below a certain level of deformation the cell behaves as an elastic material and slowly return to its initial state (Proseus et al., 1999). In case of irreversible deformation, the cell would acquire a new shape assisted by the enzymatic effect which replenishes the wall with new structural materials (Dumais et al., 2006). Moreover, plant cells are bound together by the cell wall, meaning that, contrary to animal cells, they cannot slide along each other or migrate to other regions. A consequence of this feature is that plant tissues are very rigid and single cell movements are transmitted in cascade to neighboring cells creating a tissue-spanning mechanical stress (Hejnowicz et al., 2000). As mentioned earlier, the first attempt to derive a set of mathematical equations able to describe mechano-hydraulic cell growth was proposed by Lockhart (1965), which stated that cell expansion rate is a function of cell volume, cell wall extensibility, turgor pressure, and turgor yield threshold. This equation represented the first attempt to couple water uptake and cell wall mechanics, where irreversible cell expansion is driven by the action of internal turgor pressure (Lockhart, 1965). However, the role of water fluxes is neglected in the majority of plant models, by simply assuming that the turgor pressure is a constant driving force. This unrealistic assumption is acceptable for single cell experiments but presents some practical problems in a multicellular context: neighboring cells grow at different rates  and water flow is affected by the relative geometry of the interconnected cells (Dinneny, 2020) as water availability is not uniform along the tissue (Robbins and Dinneny, 2018). Moreover, growth rate and pressure level are not always correlated even in single cells (Rojas and Huang, 2018). To reconcile Lockhart's equation with these issues and generalize it to a multicellular environment, a recent study proposed a mechanohydraulic model (Long et al., 2020), where cell growth and turgor pressure can autonomously emerge from the interaction of tissue mechanics and tissue hydraulics. Using Frontiers in Plant Science | www.frontiersin.org atomic force microscopy the authors showed significant variability in turgor pressure between cells (Long et al., 2020). Smaller cells resulting from cell division presented higher internal turgor press, while larger cells had lower pressure. In a situation where cell division does not occur, cells growth rate would tend to homogenize with turgor pressure decreasing as cell size increase (Long et al., 2020). Overall, despite its shortcomings Lockhart's equation is still today the base model for many applications in cell physiology. A major technical issue when dealing with growing organ is how to approach cell division (Figure 2A). When considering the organ as a continuous form without single-cell identities (i.e., FEM methods), cell division is typically omitted and organ growth is usually obtained through iterations of tissue remeshing and growth cycles. The problem of cell division rules had sparked interest since the nineteenth century, when Herrera proposed the so-called shortest wall solution (Errera, 1888). In 2D models, cells are commonly represented as polygons, so Herrera's rule can be implemented by simply using as division plane the shortest possible edge that divides this polygon in half (Figure 2A). The evidence seems to suggest that cell division does not occur at random (Sahlin and Jönsson, 2010), and that each cell locally controls its own division rules (Besson and Dumais, 2011). Several studies have also confirmed the existence of strong genetic control over cell division geometry (Dong et al., 2009). In plants, there is a strong correlation between the division plane and microtubule orientation (de Keijzer et al., 2014). Therefore, it has been proposed that cell divides orthogonal to the principal direction of growth where the maximal tensile stress occurs (Louveaux et al., 2016).
The combined action of cell growth and division gives rise to different cellular patterns and eventually determines the global organ shape. Even the simplest multicellular organisms may exhibit complex cellular patterning (Dupuy et al., 2010). Plants produce diverse geometric shapes, such as flowers, leaves, and networks of roots. A prominent example of patterning in plants is phyllotaxis; the process in which leaves (but also flowers and petals) are arranged around the growing vegetative stems. The beautiful geometry of phyllotaxis has attracted the attention of botanists and mathematicians since ancient times (Adler et al., 1997), and the process is known to be driven by the interaction between the phytohormone auxin and tissue growth to optimize light capture (Strauss et al., 2020). Phyllotaxis has been the subject of many attempts to formalize its mechanism through computational modeling techniques (Mitchison, 1977;Jönsson et al., 2006;Smith et al., 2006a;Zhang et al., 2021).
Morphogenesis in plants is not limited to global multicellular interplay, as whole organs are able to direct growth to form shapes designed to accomplish a specific function or to better adapt to the environment. Leaves in particular are known to exist in a plethora of shapes and dimensions, despite their almost indistinguishable primordia (Vlad et al., 2014;. Leaves emerge from the shoot apical meristem as a result of the phyllotactic process described above, and their patterning is regulated by the phytohormone auxin , which promotes growth and differentiation through the formation of maxima along the leaf margins ( Figure 2B).
Nutrients are carried along the leaf surface throughout the venation system, and the chaotic distribution of this network of veins markedly contrasts with the pristine symmetrical beauty of other processes like phyllotaxis. Nonetheless, researchers managed to devise models of venation, where the apparent complexity of these patterns is simply the result of a selforganizing process of continuous leaf surface "colonization" driven by the tandem action of auxin and cellular growth (Runions et al., 2005). Several other computational models have been developed to understand the diversification of leaf geometry (Lindenmayer, 1977;Rolland-Lagan et al., 2005;. Plants also possess radially symmetric organs such as the primary root ( Figure 2C). This rod-shaped organ has the main purpose of penetrating through soil in search of nutrients, as well as providing anchoring and stability to the plant (Petricka et al., 2012). The typical Arabidopsis root is an anisotropic structure made of radially distributed cells, with different tissues longitudinally oriented and clearly distinguishable under the microscope (Dolan et al., 1993). The root meristem in particular displays a surprisingly conserved cellular organization, with each cell type occupying a defined position and playing a determined role during cell proliferation (Van Norman, 2016). Root growth is regulated by auxin which flows through the inner tissues toward the tip, where it accumulates creating a maximum, and it is later refluxed back through the outer tissues (Kepinski and Leyser, 2005). The mechanisms behind root growth have been studied and tested in-silico using several computational models (Grieneisen et al., 2007;Fozard et al., 2013;Bassel et al., 2014;Jensen and Fozard, 2015;Marconi et al., 2021).
We have provided examples of complex plant organ shapes found in nature and in the following sections, we will present various modeling methodologies developed to address organ shape complexity. A comprehensive model of morphogenesis requires the definition of a digital structure underlying the biological tissue of interest as well as a set of biologically sound rules for growth and patterning.

DIGITAL REPRESENTATION OF PLANT TISSUES
The first problem to consider when modeling morphogenetic processes in plants is the digital representation of the underlying plant organ geometry. Whereas biological tissues are inherently continuous, the internal computer memory only allows discrete elements. The vast majority of models described in the literature typically rely on several "digital" paradigms. However, as memory and computational power are limited, tissue topology representations must balance model performance with model accuracy. The "digital" representations are concerned with a representation of plant tissue topology amenable to user interaction by providing both ease of use (sometimes with a graphical user interface), and flexibility, allowing the user to control complex biological processes such as growth and cell proliferation.

Lattice-based and Particle-based Representations
One of the main issues faced when embarking on the challenging problem of modeling multi-cellular organisms is the level of cell structure abstraction. Large tissues composed of thousands of cells can be very computationally demanding (Krul et al., 2003), and some form of simplification is often unavoidable. Lattice-based models are a derivation of the notion of cellular automata initially proposed by Ulam (1962), where cell-to-cell interactions are regulated by transition functions dependent on the current state of the cell and its neighbors. In lattice-based models each cell responds independently to external stimuli and follows cell-specific rules (Hwang et al., 2009). A common extension of latticebased models suitable for modeling tissues is the cellular potts model (CPM; Glazier and Graner, 1993). A CPM is composed of a grid-like lattice where each site is denoted as a "pixel" possessing different identities, such as the cell or the medium ( Figure 3A). Cells are allowed to grow, divide, and interact with each other, and the output of a CPM results from the interaction of the collective behavior of individual entities and not from global rules acting on the whole system (Voss-Böhme, 2012). Cellular potts models use an energy-based approach to simulate growth by minimizing the total energy of the system. A main concern with CPM is their inability to properly According to the shortest-distance rule (Errera, 1888), the cell is divided by identifying the shortest division plane cutting through the cell. To avoid conflict with existing vertices, the two ends of the division plane can be spread apart by a minimal threshold distance. Cell division splits the cell into two daughter cells yielding two new connected vertices in the cell wall. Often, the cell wall is "pinched" along the newly created vertices (Smith et al., 2006b).
(B) Computational model of leaf shape development (adapted from . Mature leaf shape is achieved through the interaction between three components: a proximal-distal hypothetical morphogen and two master regulators. The basipetal red-orange-yellow gradient region defines the action of the growth morphogen (red: higher growth, yellow: lower growth), where the dashed line indicates the border between the actively growing region and the differentiation zone. The marginal patterning of the leaf blade is the result of the combined action of a local growth activator (red) and a growth suppressor (blue).
(C) Computational model of radicle emergence (adapted from Marconi et al., 2021). This simulation reproduces the embryonic emergence of the root meristem of Arabidopsis. Organ growth follows the combined action of cell elongation, cell division, and tissue mechanics. Note that the organ maintains a distinct anisotropic form through the self-organization of cortical microtubules despite each single cell being expanded by uniform turgor pressure (see Marconi et al., 2021 for details reproduce the effect of the plant cell wall (Mele et al., 2015), but a solution to this problem was proposed with the formulation of hybrid mechanical systems (Merks et al., 2011). CPMs have been previously used to model auxin dynamics in the growing root of Arabidopsis thaliana (Grieneisen et al., 2007;Laskowski et al., 2008;Mähönen et al., 2014), leaf venation and meristem development (Wabnik et al., 2010;Merks et al., 2011), shoot apical meristem generation and maintenance (Banwarth-Kuhn et al., 2019), vascular tissue formation (Banwarth-Kuhn et al., 2019), and division plane selection in A. thaliana (Moukhtar et al., 2019). Also, non-CPM lattices have been employed to model the Arabidopsis root (Mironova et al., 2010).
Analogous to lattice-based models are particle-based models (MacAl and North, 2010), where each element (usually cells) is represented as a single particle connected to other particles by permanent bonds. The tissue system is updated by steepest descent minimization which relaxes the forces between the particles. Models of plant tissues implemented with particle systems have been mostly applied to patterning static non-growing tissues (Deinum et al., 2017;Chakrabortty et al., 2018;Mirabet et al., 2018) or agent-based modeling of auxin transport dynamics (Garnett et al., 2010), but they can also be in principle extended to simulate both growth and cell division (Van Opheusden and Molenaar, 2018). Each element (or pixel) of the grid possesses a specific identity; in this case, we have two cellular types (indicated by 1 and 2) and a medium (no labeling). CPMs are solved by minimizing the Hamiltonian energy of the system, which allows the elements of each cellular type to "group together" and isolate from the surrounding medium. (B) Two possible ways to represent multicellular tissues with vertex-based models. Cells can be identified by their cell walls (left). Cell walls are represented as polygons made of vertices connected by edges shared with adjacent cells (in accordance with the biological properties of plant cell walls). Cells can be identified by their centroid (right). Cells are represented by a network of vertices and edges connecting the centroids of adjacent cells. Cell walls and boundaries can be abstracted in several ways, i.e., using Voronoi partitioning (Mosca et al., 2018). (C) L-systems modeling (adapted from Prusinkiewicz et al., 2018). An idealized realization of an L-system simulation is shown. The different stages of development of a simple tree branching structure are obtained from simple axioms and recursively applying procedural rules to a small group of elements. (D) Schematic representation of the "cell complex." A simple two-dimensional cell complex made of two connected triangles (left) can be represented as an incident graph (right). Notice how the three cell dimensions (vertices, edges, and faces) occupy three different levels and share boundaries with each other. There are also two pseudocells ┴ and ┬ as the infimum and supremum of the incidence graph (see Prusinkiewicz and Lane, 2013

Vertex-Based Graphs
From the onset of embryonic development up until the final adult stage, plant tissues are subject to the effect of internal and external forces (Ten Hove et al., 2015). Internal forces are the result of cell expansion and proliferation, mostly driven by the action of turgor pressure (Shabala and Lew, 2002). External forces instead result from gravity and the interaction with the environment (Masson et al., 2002). Plant cells are interconnected by a common structure known as the cell wall, which means that unlike animal cells they are not allowed to freely dislodge from their current position and independent mechanical movement is highly constrained (Liepman et al., 2010). A deeper understanding of tissue morphogenesis requires taking into account the laws of mechanics acting on the cell wall. Internal and external forces induce cell wall deformation, which is difficult to model without cellular-level abstraction. Vertex models were proposed to investigate the rules governing cell motility and cell-cell interactions during morphogenesis (Weliky et al., 1991). In a vertex model, cell boundaries are represented by a network of vertices interconnected by edges (Davidson et al., 2010). Therefore, a single cell is rendered as a polygon, where each edge represents the cell wall shared with a neighboring cell or the external space ( Figure 3B).
The mechanical properties of a vertex model are usually implemented on top of the vertex-edge graph as a simple mass-spring system, where cell growth is driven by updating the resting length of the edge springs (de Boer et al., 1992). Vertex models can also include advanced formulations of the Newtonian laws of motion to create more flexible applications (Fletcher et al., 2014), and they can be easily extended to simulate 3D structures (Alt et al., 2017). Vertex models have been used as the theoretical foundation for the Vertex-Vertex (VV) system (Smith, 2006). The VV system describes a methodology for modeling dynamical surfaces on a discrete 2D manifold topology, which generalizes vertex models to make them practical in a wide variety of situations. The VV system has been released with a ready-to-use implementation in C++, and successfully used in a range of applications, such as phyllotaxis (Smith et al., 2006a;Smith and Bayer, 2009;Hartmann et al., 2019), lateral root response to gravitropism (Waidmann et al., 2019), establishment of apical-basal axis in the plant embryo (Wabnik et al., 2013), leaf shape development , leaf venation patterns (Runions et al., 2005), tissue cell polarity establishment (Abley et al., 2013), cells shape lobbiness , root growth on nutrient availability (Ötvös et al., 2021), and lateral root priming (Perianez- Rodriguez et al., 2021).

L-Systems
The need for models capable of describing elaborate morphogenetic processes in flowering plants prompted many researchers to focus on models capable of producing selfgenerating structures. Lindenmayer (1968) proposed a formal grammar system inspired by cellular automata the L-systems . L-systems are descriptive models based on natural language processing that simulate plant growth and organ development by assuming an initial set of symbols and production rules that can recursively expand the original set into more and more complex fractal-like structures ( Figure 3C). Contrary to cellular automata, which deploy a space-centered Eulerian perspective, L-systems aim at establishing a structure-focused Lagrangian view (Prusinkiewicz and Runions, 2012;Prusinkiewicz et al., 2018). The modular structure of plants therefore represents the perfect application ground for the L-system approach. L-systems are versatile parametric models that allow the incorporation of molecular-level processes and genetic regulatory networks, and they have been used to recreate the vegetative development of simple multicellular organisms like Anabaena (Lindenmayer, 1975) or complex forms such as trees (Allen et al., 2005). This technique has been successfully applied to important research topics in plant development including phyllotactic patterning (Strauss et al., 2020), epidermal cell shapes , leaf shape emergence , virtual crop generation (Marshall-Colon et al., 2017), auxin-driven patterning (Cieslak et al., 2015), control of bud activation , apical hook formation (Žádníková et al., 2016), fruit expansion , and inflorescence (Owens et al., 2016). L-systems can integrate external stimulus (i.e., temperature effect) and allow the prediction of plant phenotypes (Palubicki et al., 2009). Recent studies have also shown that L-systems can be combined with stochastic simulation algorithms to overcome typical limitations of purely deterministic model description (Cieslak and Prusinkiewicz, 2019). Furthermore, L-systems were combined with deep learning for the robust image processing of organ structures (Ubbens et al., 2018).

The Cell Complex
One feature peculiar to plants is that cells cannot move with respect to one another and the only event affecting tissue topology is cell division. The "cell complex" is a recently developed system capable of capturing the topology of plant (Prusinkiewicz and Lane, 2013). The "cell complex" is defined by mathematical elements named "cells" (not to be confused with biological cells) of different topological dimensions (i.e., 0 for vertices, 1 for edges, 2 for faces, and 3 for volumes), all of them organized into coherent structures called "flip tables. " The flip tables are sufficient to represent the whole system, and they provide all the basic working operations such as iterating, merging, splitting, getting geometric information (orientation, boundaries), and so forth ( Figure 3D). This representation is used to build the Cell Complex Framework, a C++ API which can be used for computational modeling in 2D and 3D. This API is part of the advanced 4D modeling framework MorphoDynamX. 1 The cell complex is a recent innovation and therefore its applications are still scarce; nonetheless, models using cell complexes have been applied to heterocyst formation in Anabaena and leaf margin morphogenesis (Prusinkiewicz and Lane, 2013) as well as to Frontiers in Plant Science | www.frontiersin.org cell division in the Arabidopsis embryo (Yoshida et al., 2014). At the same time, using the cell complex it was possible to recreate the complete 4D map of early Arabidopsis embryo development including all successive cell division events (Yoshida et al., 2014). Similarly, the cell complex was chosen as the base structure for the simulation of root emergence and development in a comprehensive model combining mechanical and biochemical signals .
In summary, we reviewed some popular methods used to represent cells and tissues, focusing on modeling plant development. The next step is to extend such biological representations with the specific rules applied to organ growth and biomechanics.

CELL GROWTH MODELING APPROACHES
After choosing the underlying digital representation of the plant organ one may apply a system of growth rules that allows the initial structure to evolve and generate the biological patterns found in nature. The majority of growth models focus on approximating the physical laws of movement. Organ growth can be coordinated either on global or local level; by defining general rules of organ development, or by specifying local rules for the individual elements composing the organ (i.e., individual cells). The accuracy of the representation largely depends both on the level of abstraction and the complexity of the physical rules. General biological models usually employ compact implementations that may simplify all the nuances of the physical laws of motion to provide qualitative and more universal solutions. In contrast, more specific or focused models dive into the mechanistic description of the underlying organ growth processes.

Descriptive Growth Rules
In general, it is possible to achieve complex organ shapes by specifying global regulation for the entire organ or alternatively by determining local rules at cellular level. Local growth can be difficult to define as it requires controlling mechanical interactions between local entities to avoid conflicts and breaking the laws of physics (Rodriguez et al., 1994). A common approach at solving this problem is to use descriptive models of growth (Mosca et al., 2018). This type of models is sometimes called "organ-centric system, " as growth is accomplished by transforming plant tissues along a determined path in an Eulerian fashion. Organs that display a simple conserved geometrical shape can be represented with a descriptive growth model, as cell movement can be indirectly achieved by moving the cell along a path pre-defined by a mathematical function ( Figure 4A). For example, the Arabidopsis root tip presents a radially symmetric geometrical shape resembling an elongated cylinder with a smoothed end (Dolan et al., 1993). Hejnowicz and Karczewski (1993) developed a descriptive model for root growth that is specified in linear parametric coordinates and later mapped onto the curvilinear Cartesian coordinates, this allows to easily reconstruct the root tip curvature. A similar idea has been applied to the plant shoot apex (Kucypera et al., 2017), by mapping polar coordinates to Cartesian coordinates. A major advantage of organ centric systems is that growth can be specified in only one dimension, thereby greatly reducing system complexity.

Mass-Spring Systems
In a classical mass-spring system, the edges representing the cell membrane/wall are idealized mechanical springs each connected by two masses constituted by vertices (Figure 4B; Mosca et al., 2018). Internal forces caused by the cytoskeleton are sometimes ignored as they are much weaker compared to the forces exerted on the cell wall. A damping force is commonly applied to the masses to prevent numerical oscillations around the steady state. Springs usually have negligible mass, while the mass of vertices depends on the properties of the underlying tissue. Spring can be modeled as elastic, viscoelastic, or plastic elements, depending on the nature of the biological material. For instance, tissue expansion in vertex-based models relies on internal turgor pressure, which acts on the vertices masses of the cell membrane forcing the springs to elongate, while tissue growth can be easily achieved by updating the spring resting length. Mass-spring are arguably the simplest mechanical systems to implement and have been used to model the shoot apex (Hamant et al., 2008;Uyttewaal et al., 2012), petal shape (Rolland-Lagan et al., 2003), leaf venation networks (Corson et al., 2009a), shoot apical meristem (Corson et al., 2009b), tip-growing cells (Dumais et al., 2006), shoot apex cytoskeleton (Hamant et al., 2008), and cell-cell interaction (Dupuy et al., 2008), root development (Waidmann et al., 2019), and apical hook bending (Žádníková et al., 2016). A major limitation of the mass spring system is its dependency on mesh discretization, as alternative meshes may produce dramatically different results. Moreover, complex mechanical processes such as anisotropy are challenging to implement with these systems (Mosca et al., 2018).

FEM Models
Finite element method (Becker et al., 1982) is a popular modeling technique borrowed from the engineering field and often used for modeling continuous mechanics of plant development. For many complex problems involving partial differential equations (PDE), the analytic solution is often not available, and it is necessary to recur to numerical methods to solve it. FEMs discretize a continuous mechanical problem to derive it efficiently over two or three space dimensions. Solving a FEM requires several steps that define how physical properties are applied onto a predefined geometrical structure ( Figure 4C). First, we define the initial geometry of the system, created either from segmented sample images or manually drawn using dedicated computer tools (Bassel and Smith, 2016). The geometry is later subdivided into small discrete elements by constructing a polyhedral mesh of the object of interest. Each element is defined as interconnected nodes and the global behavior of the system is constrained by boundary conditions imposed on the structure. The modeler then applies external forces and defines the interactions between the elements of the system. The mechanical rules describe how the elements behave under forces including stresses and tension (Hamant et al., 2008). Applying FEMs to model pavement cells revealed the local distribution of mechanical stresses (Sampathkumar et al., 2014). The cell wall is also another candidate system to be studied with FEMs, for example, by simulating the effect of turgor pressure over different wall shapes and thickness (Forouzesh et al., 2013). Tip growth of pollen tubes has been studied through FEM models as a hollow shell with uniform thickness (Fayant et al., 2010;Vogler et al., 2013), or to quantify data from micro-indentation (Bolduc et al., 2006). Similarly, FEMs were employed to simulate the branching morphogenesis of Arabidopsis trichomes (Yanagisawa et al., 2015), revealing a  Smith et al., 2006b). The shoot surface layer can be generated by a B-spline rotating around the longitudinal axis of the shoot apex (upper right plot). A point p located on the surface of the apex is located at coordinates (θ, a) moves away from the apex tip with the velocity v(a), where a is the distance from the apex tip, measured along the curve on the apex surface. The shape of the generating curve and the angle θ determines the growth in Cartesian coordinates. The two sample plots show the relative elemental rate of growth (RERG) as functions of distance from the apex tip, either at constant growth (lower right box, upper plot) or with a decreased growth rate at the apex tip and flank (lower right box, lower plot). (B) Schematic of a cell represented by the mass-spring system. The cell wall is a polygon whose vertices and edges are masses and springs, respectively. The effect of internal turgor pressure generates forces perpendicular to the wall edges, and distributed over the vertex masses. The vertices movement is further restricted by the elastic forces exerted by the springs. (C) Schematic representation of the finite element method (FEM). A continuous object for which a global function of growth needs to be approximated is subdivided into a conglomerate of smaller geometrically manageable finite elements (triangles in this case) in a process called "meshing." Each element is represented by a set of piecewise linear equations (sometimes called "trial" functions) derived from the original problem. After the calculations are done for each finite element all sets of equations are systematically recombined into a global system that approximates the original one. Generally, a finer meshing results in a better approximation. (D) Comparison between the force-based approach and position based dynamics (PBD). In a force based approaches the total cumulative force is calculated for each physical entity. For instance, a mass-spring system includes the total force acting on the vertex via internal turgor pressure and the elastic forces of the springs. Acceleration and velocity are calculated for each time step (usually recurring to explicit or implicit integrating solvers), and the object position is finally updated. In the PBD approach, the object is subjected to different constraints instead (in this case, two distance constraints and an area constraint, see Müller et al., 2007;Marconi et al., 2021), and each constraint is sequentially projected on the object updating its position. The vertex velocity is recalculated in the final step.
Frontiers in Plant Science | www.frontiersin.org strong axial growth caused by the transversal alignment of microtubules. In the context of shape deformation, FEMs were applied to reconstruct the reversible shape changes of stomata guard cells of leaves pores (Bidhendi and Geitmann, 2018), to describe the mechanical feedback restricting sepal growth and shape (Hervieux et al., 2016), and the emergence of epidermal cell shapes . Robinson and Kuhlemeier (2018) constructed a finite element model to understand the effect of external stresses on the hypocotyl, while a study of Bassel et al. (2014) recapitulated the mechanical outgrowth of an Arabidopsis root radicle. At a higher structural level, FEMs have been used to explain the modality of explosive seed dispersal (Hofhuis et al., 2016). A powerful feature of FEM modeling is the integration of quantitative experimental data. Live cell imaging of the shoot apical meristem combined with finite-element modeling allowed to unravel the functional distinction between differential growing tissues beyond the sheer genetic specification (Kierzkowski et al., 2012). FEMs can also be combined with other modeling techniques; Fozard et al. (2013) (Jensen and Fozard, 2015) proposed a hybrid mass-spring/FEM "vertex-element" model by merging a vertex-based structure with a finite-element discretization of in-plane walls; while Boudon et al. (2015) succeed in applying FEMs to mechanical models of young meristems at "pseudo"-cellular resolution. Finally, understanding the relationship between gene activity and organ shaping led to the development of third-party computer applications such as GFtbox (Kennaway et al., 2011), which aims at integrating different components of plant morphogenesis into a coherent software interface based on FEM.
Overall, finite element modeling is a powerful tool that allows researchers to simulate the behavior of complex plant organ shapes and their interactions with genetic activity and environmental stresses; however, at the cost of high computational demand (as the organ grows) and (usually) lack of singlecell definition.

Position Based Dynamics
Force-based systems (such as FEM and mass-spring models) represent the typical approach to deal with the mechanical properties of plant systems. Under this paradigm and Newton's second law of motion, acceleration is computed from total internal and external forces. An explicit or implicit time integration method is then used to update the velocities and the positions of the mesh elements. Unfortunately, force-based systems are often unstable and require slow implicit integration to converge, which is a major limitation for modeling complex organ structures at cellular resolution.
Physically-based animation has always been a major topic in computer graphics (Bargteil et al., 2020). This field is concerned with finding new methods for the simulation of physical phenomena such as rigid body dynamics, objects deformation or fluid flow, with a strong focus on stability, speed, and robustness. Position based dynamics (PBD) is a recent physically-based animation technique that prevents the typical instability problems of force-based systems through the action of local constraints (Müller et al., 2007). With this approach it is possible to omit the velocity layer and immediately work on positions, making the system more stable and controllable ( Figure 4D). PBD has been applied in a large number of applications, mainly outside the biological world (Bender et al., 2017). Our own study  has shown its potential benefits in plant morphogenesis by implementing PBD inside the "cell complex" environment (Prusinkiewicz and Lane, 2013) and using it to model the growth or the Arabidopsis root . Similar to other modeling techniques developed for physically-based animation, PBD provides higher performance and stability, thanks to the aforementioned replacement of force-based motion with ad-hoc designed constraint functions. However, as expected, such advantages come at the price; Newtonian forces are abstracted through constraints and therefore the dynamics of the system do not have a direct physical interpretation (Müller et al., 2007). Another limitation of PBD is its dependency on time step size and number of iterations (Bender et al., 2017). To address these problems, an alternative formulation of PBD has been derived by introducing a new constraint formulation that corresponds to a well-defined concept of elastic potential energy, which allows for solving constraints in a time step and iteration count independent manner (Macklin et al., 2016).
Taken together, modeling of organ growth is a complex endeavor as it requires avoiding conflicts between moving elements such as cells or tissues, while preserving the global structure of the organ. However, combining growth mechanics with biochemical reactions represents an increased level of complexity and is currently an active research subject.

SURVEY OF MECHANICAL AND BIOCHEMICAL MODELING FRAMEWORKS
The applications of computational techniques in biology have fueled the development of extendable frameworks to avoid the multiplication of models tailored for just a single problem, providing a fully integrated environment for computational modeling and hypothesis testing . In many instances, these frameworks can provide the user with a fully operational toolkit that abstracts the underlying implementation and allow the user to focus the effort on addressing the actual biological question. Several software packages have been developed to satisfy different needs for plant-related scientific problems, and most of them offer open-source licenses.
L-studio is a Microsoft Windows software that uses L-system (discussed above) to simulate models of plant morphogenesis; the corresponding version for Linux is known as Vlab (Karwowski and Prusinkiewicz, 2004). This modeling package represents one of the oldest software specifically designed for plant structure modeling. L-studio objects are C++ modules that can be loaded and executed into the main simulator. The simulator produces a visual representation of the model using the OpenGL graphics library. 2 L-studio also provides a browser for organizing and accessing objects such as plant segments, cells, or tree branches on local and remote machines as well as a series of editors and other modeling tools for creating and modifying these objects.
OpenAlea is an user-friendly open-source platform that provides the researcher with a graphical user interface comprising a set of tools specifically dedicated to plant tissue modeling (Pradal et al., 2008). OpenAlea was designed to be an easy to use, reusable and extendable collaborative environment. The main advantage of this software is the definition of the model using a graphical language. Model components can be visually edited by the user and connected to other components. Each component contains rules and parameters that biologically define the execution of the model. Users can add new functionalities to the system as Python 3 scripts through the package manager at run-time without modification of the framework. The current package includes a number of readyto-use external modules like the VPlants package for plant architecture analysis and the PlantGL graphic library for plant geometric modeling. OpenAlea has been used in a number of models to describe processes ranging from auxin transport to root branching (Lucas et al., 2008;Perrine-Walker et al., 2010;Rutzinger et al., 2011).
VirtualLeaf is a cell-based computer modeling framework designed for plant tissue morphogenesis (Merks et al., 2011). This program allows to model cells, cell walls, chemicals diffusion, and reactions, and to define rules that regulate growth and development. VirtualLeaf inner mechanism is similar to the previously described CPM. One of the main motivations that justified the creation of this software library was the observation that plant cells are not allowed to slide along each other but are instead constrained by the presence of the cell wall. Common CPMs cannot prevent this issue, while VirtualLeaf offers an alternative off-lattice method that implements a working solution. VirtualLeaf has been successfully used to aid understanding plant morphogenesis in several different contexts (Wabnik et al., 2010;van Mourik et al., 2012;De Rybel et al., 2014;Qi et al., 2017).
Multicellular organisms are capable of producing extremely different shapes coordinated by interaction between individual cells. Cells communicate by different means through the exchange of mechanical and chemical signals. The desire to explicitly represent this network of interactions has driven the creation of CellModeller (Dupuy et al., 2008). CellModeller is a generic software tool designed for the analysis and modeling of plant morphogenesis at cellular resolution. This framework can execute systems with more than 1,000 cells and their interactions. Moreover, it is specifically designed for plant tissues, and it can be easily expanded with additional models using the XML file format (Dupuy et al., 2008).
The popular finite element method (described in the previous section) is the mathematical background for the GPT-framework (Kennaway et al., 2011). GFtbox, is a MATLAB application devised for flat organs commonly such as leaf and petal, but it is also compatible with other organ forms. Different patterns of growth can influence the deformation of continuous sheets of biological tissues creating vastly dissimilar shapes. The GFtbox can combine mechanical growth with gene activity, tissue identity, and biochemical properties at the organ level. For example, the GFtbox has been successfully implemented to describe growing polarized tissues and genetic control of morphogenetic processes (Kuchen et al., 2012;Sauret-Güeto et al., 2013).
MorphoGraphX is one of the most popular tools for visualization and analysis of 4D biological datasets (de Reuille et al., 2015). It allows the user to extract organ/cell shape through segmentation and to quantify cell growth and signal fluorescence. Moreover, this software allows the dynamic modeling of templates extracted from imaging data through simple scripting. MorphoGraphX can be easily expanded with external plugins. One of them is MorphoRobotX, an add-on used to analyze data from the cellular force microscope, a custom microindentation system that is specialized for plants . Recently, the new-generation version of MorphoGraphX called MorphoDynamX has been constructed based on the "cell complex" data structure described in the previous section (Prusinkiewicz and Lane, 2013). This addition greatly improves the opportunities of modeling subdividing geometry in 2 and 3 dimensions, a big step forward in comparison with simpler vertex-to-vertex technologies. MorphoDynamX and the cell complex are very recent developments but they have already shown their potential in the modeling of cell division in the Arabidopsis embryo (Yoshida et al., 2014) and the organogenesis of Arabidopsis root . MorphoDynamX can also be expanded with external modules; for example, MorphoMechanX is an add-on for MorphoDynamX that enables the mechanical modeling of biological tissues using FEMs for solids elements (Long et al., 2020;Natonik-Białoń et al., 2020;Läubli et al., 2021).
Outside of plant world other frameworks could be adapted to model plant organs such as CompuCell3D (Cickovski et al., 2005) -a lattice-based multicellular modeling environment to simulate tissue formation by merging the CPM with chemical fields and diffusion equations. CompuCell3D functionalities have been demonstrated in the context of animal embryogenesis, cell population dynamics, tumor formation, and many more (Andasari et al., 2012;Swat et al., 2012;Fortuna et al., 2020;Liu et al., 2021). Another valuable tool to simulate animal embryogenesis is MecaGen, a C++ simulation platform of animal multicellular development relying on mechanistic agentbased models (Delile et al., 2017). This software is capable of combining cell mechanics with gene expression and intracellular signaling. For example, this tool has been applied to replicate zebrafish epiboly collective cell behavior (Delile et al., 2017), epithelialization of Drosophila wings (Neto-Silva et al., 2009), and cell fate specification in animal tissues (Chan et al., 2017).
To date, diverse computer modeling frameworks for plant organogenesis have been developed that typically integrate either biochemical or mechanical cues. These frameworks were critical in providing answers for many important scientific questions that otherwise could not be addressed with purely experimental Frontiers in Plant Science | www.frontiersin.org approaches. However, there is still a shortage of user-friendly model frameworks that combine mechanics with gene regulation and intracellular signaling at single cell level, thus leaving a vast space for future improvements.

CONCLUSION AND PERSPECTIVES
Nowadays, computer models are a fundamental tool in virtually any field of scientific research. In plant science, the quantitative representation of organs and tissues still remains computationally challenging, in particular given the amount and quality of data needed to resolve growth at the cellular level. Moreover, constructing advanced models of morphogenesis requires the researchers to be competent not only in experimental techniques but also to possess a basic knowledge of programming, computer graphics, physics, mechanics, and biostatistics. Still, the potential benefits of including computer models into experimental research environments are vast; many are the cases where naive intuition about a specific natural process falls short, since deterministic consequences and causal connections that appear logically obvious can in fact be dangerously misleading. Computer simulations have demonstrated that complex forms can emerge from the dynamic interactions between minimal components of the system. Therefore, models allow researchers to separate the inquiry from his cognitive bias and to test hypotheses in an objective and reproducible way.
However, we must be extremely careful when interpreting and reporting the results obtained from computer simulations. Just like any other analytical approach, models are subject to the evergreen motto "garbage in, garbage out"; by feeding the model with erroneous, biased or imprecise data, we can at most expect the same quality as output. There are other common mistakes in computer modeling, related to the model definition: extremely simple/complex models that tend to underfit/overfit the data, low-definition graphical representations of the biological structure, unrealistic physical or chemical properties regulating system development, lack of practical use and self-fulfilling expectations fueled by bias reinforcement, among others. A strict collaboration between computer scientists and wet-lab researchers is of paramount importance to guarantee the biological relevance of computer simulations. Additionally, a continuous feedback between the model and experimental data ensures that the necessary assumptions and simplifications built into the model do not compromise the final inference.
In the area of development biology, these complications are further exacerbated by the intrinsic sequential nature of the subject under scrutiny: growth occurs over time and often specific developmental stages are hard to observe or completely unknown, this produces holes in the idealized mechanistic process that need to be carefully filled in accordance with biologically sound principles. Recognizing the role of spacetime interactions and the importance of biomechanics in plant growth is fundamental to achieve a realistic simulation, as modeling can be expressed in 1D, 2D, and 3D (+ time). The modelers must also decide which level of abstraction should be considered; this means, for example, choosing whether to represent the tissue as a continuous material rather than subdividing it into single cells or other type of discrete representation. Depending on the process under investigation either approach can be more appropriate than the other.
The recent years have witnessed the development of a plethora of computational tools designed to unravel a disparate range of problems in plant morphogenesis. Nonetheless, there still remain many unanswered questions; what kind of feedback exists between cell division and mechanical/molecular processes? How do cells regulate their own anisotropic growth? What kinds of forces determine cell polarity (e.g., how are hormones transported and how are they distributed in a polar manner)? How does the cell wall expand under internal turgor pressure but withstand external forces? How do organs achieve a desired shape without the aid of global long-range signals? These are just a few open questions of which we possess only a partial knowledge, many others laid unanswered. We have the feeling that the future of computational modeling is bright and these models will yet play an essential important role in deciphering the mysteries of plant morphogenesis.

AUTHOR CONTRIBUTIONS
MM and KW wrote the entire manuscript with no additional contribution. All authors contributed to the article and approved the submitted version.