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, 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 (Kierzkowski and Routier-Kierzkowska, 2019) 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.

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 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 (Hong et al., 2018) 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 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; Kierzkowski et al., 2019). 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 (Bayer et al., 2009), 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 self-organizing 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; Kierzkowski et al., 2019).

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.

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.

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.

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 (Routier-Kierzkowska and Runions, 2018). 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.² 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³ scripts through the package manager at run-time without modification of the framework. The current package includes a number of ready-to-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 (Routier-Kierzkowska et al., 2012). 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 (Marconi et al., 2021). 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 agent-based 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 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.

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 space-time 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.

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

This work was supported by the Programa de Atraccion de Talento 2017 (Comunidad de Madrid, 2017-T1/BIO-5654 to KW), Severo Ochoa (SO) Programme for Centres of Excellence in R&D from the Agencia Estatal de Investigacion of Spain [grant SEV-2016-0672 (2017–2021) to KW via the CBGP]. In the frame of SEV-2016-0672 funding MM is supported with a postdoctoral contract. KW was supported by Programa Estatal de Generacion del Conocimiento y Fortalecimiento Cientıfico y Tecnologico del Sistema de I+D+I 2019 (PGC2018-093387-A407 I00) from MICIU (to KW).

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.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.