Single-cell approaches for understanding morphogenesis using Computational Morphodynamics

In multicellular organisms cells grow, divide and adopt different fates, resulting in tissues and organs with specific functions. In recent years, a number of studies have brought quantitative knowledge about how these processes are orchestrated, shedding new light on cells as active and central players in morphogenesis. We explore recent advances in understanding plant morphogenesis from a quantitative perspective, defining the research field of Computational Morphodynamics. The focus is on studies combining theoretical and experimental approaches integrating hypotheses of how molecular and mechanical regulation at the cellular level lead to tissue behaviour. Finally, we discuss some of the main challenges for future work.


Introduction
Morphogenesis, the spatial development of tissues with specialized structure and function from populations of undifferentiated cells, is one of the most fascinating and complex problems in nature.Formation of functional tissues requires a strict spatiotemporal control of cellular morphology and gene expression.To achieve this, individual cells can read a multitude of microenvironmental cues, respond via their gene regulatory networks, and ultimately, undertake the appropriate fate decisions and morphological changes [1].A detailed characterization of the factors influencing cellular fate decision, the mechanisms by which they are executed at the cellular level, and how they are coordinated in multicellular tissues, is paramount to understanding the dynamic pattern formation driving development.
Single cell studies in unicellular organisms and cell cultures have made great progress in providing a quantitative description of the signalling and gene regulatory mechanisms that underlie cell fate decisions [2][3][4][5].These experiments typically involve timelapse microscopy methods, in which the cells of interest can be tracked as they grow in controlled microenvironments.Expression dynamics of key regulators can be characterised by computational quantification of fluorescent reporters [6,7] and correlated to cellular variables such as size, division and death [8].
Progress in the development of reporters, microscopy technologies and computational methods now allows for expanding these approaches to the study of developing tissues in multicellular organisms in vivo.In plant biology in particular, morphogenesis has been the object of multiple studies over the years and the role of short-range biochemical signals, long-range hormone transport, and mechanical forces in shaping tissue formation is well documented [9][10][11].The ability to quantify the spatiotemporal dynamics of gene expression and cellular behaviour at high resolution in single cells enables a finer quantitative appreciation of how signalling factors ultimately lead to specialised three-dimensional patterns.This level of detailed analysis requires powerful experimental and computational tools in order to identify, track and quantify variables of individual cells in a three-dimensional tissue as it grows in vivo.As these tools become increasingly available, it has been possible to combine high-resolution time-lapse microscopy, image processing, quantification tools and computational models in multidisciplinary efforts to better understand plant development.This integrative approach is the cornerstone of the emerging field of Computational Morphodynamics [12][13][14] (or Systems Morphodynamics [15]), which, by iterating between experimental design, quantitative data analysis and computer simulations, aims to understand the factors that shape tissues and how they affect individual cells (Fig. 1).
In this chapter, we will provide an overview of the most recent developments in the Computational Morphodynamics field in plants, particularly in terms of image processing, quantitative data analysis, and computational modelling techniques.We will review recent studies that applied these methods and technologies, with strong emphasis on work performed at the single-cell level in plant systems, taking advantage of plantspecific properties such as lack of cell migration and relatively slow growth.To illustrate some of the techniques, we will refer to several examples using Arabidopsis thaliana as biological model, particularly where single cell descriptions lead to behaviour at the scale of small tissues.We will also discuss the main hurdles and challenges in the field, as well as new methods and technologies that open new and exciting avenues for further understanding the dynamical principles and mechanisms underlying tissue formation in plants.

Capturing single cell dynamics in space and time
The accurate characterization of how single cells are affected by microenvironmental signals and respond in accordance during plant development requires a detailed spatiotemporal description of the tissue as it grows.This description is typically based on time-lapse confocal microscopy of plant tissue in vivo, or live imaging, which allows the non-destructive sampling of the tissue in three-dimensions over time.This approach, particularly in combination with other cell and molecular biology tools such as fluorescently tagged reporters, can allow for the quantification of a wealth of variables related to the dynamics of tissue mechanics, cell geometry and cell division, as well as gene expression.Such experiments require high resolution imaging data both spatially and temporally.Ideally, the number of confocal slices spanning the tissue should be as high as possible (thus covering as small a distance as possible) for individual time points, with time intervals being kept to a period as small as possible.A limiting factor in both cases is the phototoxicity that arrives from prolonged and repeated exposure of the tissue to the laser, and so a suitable compromise needs to be found.Sampling times may also need to be adjusted depending on tissue growth rates, such that individual cells can be unequivocally tracked.
The questions of how and when cells divide in plant tissues have attracted much attention.The first studies discussing cell division rules date back to the late 1800s, in which it was proposed that the geometry of the cell would determine the orientation of the cell division plane (more details on this topic in [24,25]).In the last decades, improvements of microscopy technologies connected to fluorescent markers have enabled revisiting these questions in more detail.In a pioneering study, Laufs et al (1998) used propidium iodide to stain DNA and visualise individual nuclei in Arabidopsis thaliana inflorescence shoot apical meristems (SAMs) [16].This method allowed for a detailed morphometric analysis that included the quantification of size, spatial distribution and mitotic index of individual cells in wild type and mutant lines.Propidium iodide staining effectively arrests growth, making it impossible to collect dynamical data.Reddy at al (2004) circumvented this issue by using fluorescent reporters for components of the plasma membrane, histone markers and mitotic cyclins, thus being able to observe the dynamics of cell and nuclear division events, as well as to capture cells about to or in the process of division [17].
As imaging technologies and computational power increase, higher resolution data can be generated which affords automated identification and tracking of individual cells, as well as higher precision in the quantification of variables of interest [18][19][20][21][22][23].The implementation of a powerful automated 4D imaging pipeline allowed Willis et al to observe in great detail the dynamics of cell growth and division in Arabidopsis SAMs [18].By sampling meristematic growth every four hours over the course of 2-3 days, and tracing each individual epidermal cell in space and time, the authors could show that cell size regulation does not strictly follow the so-called sizer nor adder models that had been previously described [26], but instead is an intermediate between these two paradigms.Jones et al (2017) also approached this question with single-cell quantification of cell division coupled with a mechanistic model of cyclin-dependent kinase (CDK) activity [20].The study pointed to the existence of regulatory mechanisms that dynamically regulate cell size as a function of a number of factors including growth rates.Both these studies suggest that cell size in the meristem is not an intrinsically defined (or measured) property.Beyond the timing of division, it is also of interest to understand how cells define the orientation of division planes and localization of new cell walls.Factors such as cell size, geometry and mechanical forces have all been deemed of potential relevance over the years.Recent quantitative studies have supported the view that cells divide along local minima of plane area leading to equal-sized daughter cells [21], and the alternative view that new division planes establish along directions of maximal tensile stress [22].Shapiro et al (2015) implemented a quantitative model that predicts the localization of new division planes at minima of a potential function that incorporates these and other division rules [23].Despite all its potential, quantitative image analysis at the single cell level is experimentally and computationally demanding and requires careful planning in order to avoid misleading conclusions due to technical artefacts at all stages, from image capture to quantification and analysis [27].Continuous collection of 4D data in vivo is extremely challenging in plants, due to plant growth and movement, particularly in tissues such as the SAM where plants have to be kept in light-and temperaturecontrolled chambers between sampling times, and where the imaging is usually done with water-dipping lenses [17,18].These challenges can in themselves be the source of artefacts, and the development of imaging technologies that circumvent them by allowing automated image capture with minimal disturbance can be of paramount importance.To this end, von Wangenheim et al (2017) developed a confocal microscope setup that allows vertical imaging of plant root tip growth with automated adjustment of sampling positions between time points, facilitating the tracking of objects of interest [28].This platform allowed continuous observation of root growth even when inducing rapid changes in gravity or light, and was successfully tested in zebra fish embryos, showing its general applicability.
The integration of microfluidic technologies for tissue live imaging holds great potential and, indeed, has already proved of importance in different tissues [29].One prominent example is the RootChip [30], which allows continuous imaging of different growing plant tissues while being exposed to different environmental perturbations [30][31][32].This system has been used for studying for instance gibberellin response in growing hypocotyls [31] and root-bacteria interactions [32].Microfluidic technologies afford working in very small volumes as well as the minute control over cellular microenvironment, and hence they have also become an attractive system to observe and manipulate plant cell cultures, either maintained in stable cultures [33] or freshly extracted from the tissue by cell wall digestion [34].Recently, microfluidic chips have been used to explore mechanisms of regulation of cellular geometry by applying directional auxin gradients to single BY-2 tobacco cells [35].In the same cell system, optical tweezers had previously been used to create cytoplasmic protrusions in single cells, thus allowing the charac-terization of actin and myosin dynamics and their impact on cytoplasm stiffness [36].
For some tissues, it is technically very challenging to perform single-cell quantification while maintaining the plant alive, in which case fixing tissues can be a suitable option [37][38][39].Although this makes it impossible to collect dynamical data, it allows for detailed imaging, and statistical methods can then be applied to extract information on growth and division properties from these well-defined cell patterns (Fig. 4C-D) [37].By design, confocal microscopy limits the size of the biological tissue that can be imaged at high resolution.In that regard, whole tissue and even whole organism techniques have been deployed to quantitatively characterize 3D large-scale phenomena in plant development [40,41].
While the challenges of imaging multicellular and slow growing tissues are many, several of these challenges have recently been overcome and plant development is now approaching a stage where it can be studied with live imaging tools similar to those used for bacteria and yeast, both in vivo and in cell cultures.

Quantitative image analysis
Once image acquisition has been achieved, ideally resulting in a high spatiotemporal resolution time lapse of tissue growth, quantification of cellular variables of interest is performed by deploying a set of computational methods that broadly allow the identification and tracking of each individual cell in space and time.It should be stressed, however, that image capture and analysis are not independent steps and that particular requirements of image processing should be taken into account when planning time-lapse experiments.Microscopy parameters such as laser intensity, voxel size (in particular, defining the thickness of each confocal slice), and resolution can be optimised by iterating with the image analysis software in single images before per-forming the full time course.Time sampling should also be carefully defined such that cell tracking can be performed with maximal efficiency.A number of preprocessing steps can also be crucial to improve image quality for the purpose of quantification.Deconvolution may be of importance when processing 3D images, in order to circumvent artificial stretching in the Z-direction caused by the point-spread function [42,43].Another relevant factor is plant growth that occurs during image acquisition outside the region of interest.In this case, a fast confocal scan with few slices covering the whole tissue may be acquired first, and used to correct potential growth effects in the longer acquisition, by correcting z-direction slice thickness [18].In addition to this, other operations may be beneficial such as denoising, edge detection or binarisation of the images to facilitate computational identification of regions of interest (Fig. 2A).
The identification of individual cells in a tissue is achieved by applying computational algorithms that make use of specific features of the imaging data to delimitate each region of interest (e.g. a cell, nucleus or other cellular subcompartment) in a process called segmentation.One of the most popular algorithms for cell segmentation is the watershed algorithm [44].This is particularly true for cases where fluorescence intensities are maximal at the cell membrane, with the inside of each cell having intensity close to zero.Different variants of the watershed algorithm have been successfully applied to the identification of single cells in plant tissues, from embryos [37], to SAMs [45,46] (Fig. 2B), and root meristems [45].When the fluorescent signal accumulates in a subcellular compartment such as the nucleus, other algorithms may be more suitable.An implementation of the steepest gradient ascent in particular has been used to segment individual meristem cells [47,48], and nuclei in time lapses of sepal growth [49] (Fig. 2A).Once a region of interest has been segmented, its size and morphological parameters can be quantified.The same is true for fluorescence intensity within each region, which is normally used as a proxy for mRNA or protein levels for each cell (Fig. 2C).Although this is often not the case, ideal experimental design should include a second fluores-cent reporter for a constitutive gene, which can be used to normalize fluorescence intensities in individual cells.
Once individual cells are identified, characterised and computationally catalogued for each time point, temporal correspondence must be resolved by tracking each cell, as well as their progeny, through the time course.This process thus needs to consider cell divisions as well as tissue growth and orientation changes between time points.Due to the existence of rigid cell walls, cell movement is not as much of a factor in plants as it can be in animal tissues.Optimal cell-cell pairing between consecutive time points can be achieved through a process of registration, where one of the images is linearly or non-linearly transformed to maximise matching with the other image of the pair [50].Different registration methods can be applied with reasonable success, particularly when taking into account the specificities of the tissue and objects of interest.These algorithms are available through implementations in different image analysis packages, often with convenient graphical user interfaces and the possibility for manual correction of segmentation and tracking errors [46,[51][52][53][54].In plant tissues specifically, block matching registration algorithms [55][56][57] have recently been used to successfully perform tracking of individual cells in growing SAMs (Fig. 2D) [18], and individual nuclei in growing sepals [49].
The current multitude of image processing tools provide a large pool of options that in principle should cover the exploration of a wide range of biological questions.On the other hand, such diversity also brings about the fragmentation of the community regarding tools of choice and, more importantly, file and data formats.Advances in the state-of-the-art allow us to observe biological phenomena at different scales and complexities from single molecule to whole organs [58,59].It is of paramount importance to agree on file and data standards within the community that maximise transferability of both raw image and processed data, such that we can increase power, accuracy and efficiency of quantification protocols.
The vast volumes of spatiotemporal single cell data collected through quantification of time-lapse microscopy experiments can be explored by computational methods in a considerable number of ways.An increasingly used approach is the application of machine learning methods to extract interesting features or behaviours [60].In supervised machine learning, objects of interest (e.g.cells) can be manually classified within a set of classes a priori (e.g.cell types) given a set of variables (e.g.size or morphological parameters).This dataset is then used to train a model that, once applied to a new dataset, can be able to predict with high accuracy the class of a given object by taking into account its set of variables.There are a number of models that can be used for classification such as regression methods, support vector machines, artificial neural networks, decision trees or random forests [60,61].In all these methods, training involves minimizing a function that quantifies the error between prediction and known outcome.A good model should be able to learn from the known data structure and successfully generalise to large sets of new data without overfitting.During the image processing stages, classification methods that take into account cellular morphology and tissue growth can be used to automatically perform and improve segmentation and tracking [61,62].Cellular size and morphology parameters have indeed been used for classifying cells in a number of ways, from cell cycle stage [63] to cancer activity [64].Gene expression data are also a fertile ground for the application of classification methods in order to gain biological insight.Single cell gene expression differences between related cellular populations have been used to predict the identity of genes [65] or the expression thresholds for a specific gene [66] likely to be involved in the decision to choose a particular cellular fate.In plants, a similar approach was recently followed to infer a mechanism by which expression maxima of the transcription factor ATML1 above a certain threshold during the G2 but not the G1 stage of the cell cycle predict with high accuracy the giant cell fate in growing sepals (Fig. 2E) [49].Unsupervised machine learning methods allow exploring features of interest in the data structure itself without having to a priori determine the different classes.These include clustering as well as dimensionality reduction methods such as principal component analysis or multidimensional scaling [60].
Beyond gene expression and single cell morphologies, it is of great interest to quantify spatial cellular patterns within a specialised tissue, as a means to further understand its morphodynamics.Methods that allow quantification of tissue order as a function of the spatial distribution of the cells that compose it can be of interest to understand its development and function [67].In plants, this approach has been followed to understand pattern formation in leafs, in particular the variability of trichome patterns on its surface [68].More generally, a recent study provided a topological characterization of complex plant organs by applying quantitative network analysis methods to high resolution descriptions of cellular interactions (Fig. 2F) [69].As the volume of data accumulates, it becomes crucial to integrate these measurements in order to understand how spatial cellular patterns correlate with the underlying gene expression and cellular division patterns during morphogenesis [39].
Great progress has recently been made in the ability to segment, track and quantify single cell information from 3D live imaging data.The continuous development of imaging technologies will predictably lead to data with higher spatiotemporal resolution gathered at quicker rates.The expected magnitude of data produced by these technologies will benefit from current and improved computational methods, but will also require more automatised (and preferably standardised) analysis and storage protocols in order not to become a major bottleneck in future studies.

Modelling tissue morphogenesis from the bottomup
Tissue patterning results from the interplay of intracellular regulatory networks, cell-to-cell interactions, cell growth and division, and mechanical forces [1].How can mathematical and computational models take these different processes into ac-count?And how can models integrate and be compared with experimental data?Regulatory networks have mostly been modelled using deterministic rate equations of cellular concentrations of key components.In this context, it can be assumed that, for a given set of model parameter values, the dynamics of regulatory networks will achieve certain steady states, which can be associated to different cell states or fates [47,[70][71][72][73][74].For example, different regulatory networks have been modelled for understanding epidermal patterning of hair and non-hair fates in both roots and leaves [70,71].Also, a model for brassinosteroids signalling in the root has proposed two alternative states that would represent a quiescent and a division state of the quiescent centre (QC) cells [73].
In different plant developing tissues, stochasticity in the key regulators has been proposed to play a role in cell fate decisions [75,76].In multicellular models, stochasticity has mostly been implemented as small initial cell-to-cell differences in the signalling regulators, which are deterministically propagated through the regulatory network.Though, due to the probabilistic nature of chemical reactions within cells, certain regulatory networks can show stochastic effects throughout time, and hence, stochasticity needs to be modelled in a dynamic manner.The Gillespie algorithm is an exact and discrete method that can be used for modelling the stochastic nature of regulatory networks [77,78].Yet, in the context of multicellular modelling, discrete methods that simulate the number of molecules can become computationally very expensive, and the chemical Langevin equation, which is a continuous concentration-based approach, seems more applicable [79,80] (Fig. 3).Although not being exact, this algorithm can give quite good approximations to the exact Gillespie method when the number of molecules is not too low [79].The chemical Langevin equation contains a deterministic contribution and a stochastic contribution, and the stochastic contribution becomes smaller when the number of molecules is high [79].This model formulation facilitates the possibility to use dif-ferent analytical tools that are used in deterministic systems, as nullcline and bifurcation analyses [81].
Throughout a developmental process, both deterministic and stochastic processes may simultaneously take place.Hence, to model this, hybrid algorithms are necessary, in which both stochastic and deterministic variables are integrated.In the context of giant cell fate commitment in the developing sepal, ATML1 fluctuations were simulated using chemical Langevin equations for ATML1 itself and its downstream target together with a timer variable, while cell growth was modelled deterministically [49].
Often, the patterning process happens while cells are growing and dividing.In the root, the cellular growth and division patterns seem deterministically controlled, leading to stereotypical files of cells that robustly conform the root (Fig. 2C) [73,82].The root shows a distinct behaviour; growth is mainly anisotropic, cell divisions are most importantly happening in a single dimension along the length of the root.In contrast, the SAM epidermis undergoes isotropic growth and cell divisions are occurring in multiple, but anticlinal, directions, preserving the 2D layer (Fig. 2D) [18].In the case of the developing sepal, a pattern of giant cells gradually forms while cells anisotropically grow and divide [49].Also, in leaves, the stomata pattern emerges as a result of an interplay of growth, cell divisions and cell fate decisions [83].Therefore, incorporating cell growth and cell division in the gene regulatory network models becomes fundamental for reproducing the different behaviours and spatial organisations found during morphogenesis of different tissues.
How can we incorporate growth into the models?Cell growth is the result of the internal turgor pressure and mechanical cues [9,84].In turn, signalling -and hence, regulatory networks -can modify cell wall properties, which feed back into how cells respond to mechanical cues.Turgor pressure can be effectively included into models in different ways; for instance, by imposing exponential tissue growth [49] (Fig. 3), or it can be included in the equations accounting for mechanical cues [84,85].In simulations, mechanical cues can be effectively included between vertices through spring interactions acting according to effective potentials [24] or through energies that need to be minimised [86].Alternatively, stress-strain relations (Fig. 3) can be included in triangulated cell walls, using finite element method based approaches in 3D [87].The inclusion of mechanical cues in the models and interactions between joint walls will lead to different growth rates for different cells.Indeed, heterogeneity in cellular growth rates due to mechanical factors and in asymmetric cell divisions have recently been reported [18,86].Incorporation of cell growth in the models should drive dilution effects in the different modelled concentrations, and distribution of individual molecules between daughters at cell division.If cell growth is exponential, dilution can be included as effective degradation rates of the modelled concentrations (Fig. 3).Though, given the heterogeneity and time-dependent cellular growth rates due to mechanical cues, it becomes more convenient to compute dilution effects given the actual volumetric change [49] (Fig. 3).
Taking into account the multitude of cell division patterns found in plants, another challenge is to include cell division in plant tissues models.A recent study showed cells in the SAM seemingly divide following a rule that is in between critical size and critical increment paradigms [18].However, in many cases, assumptions and simplifications of the division process will need to be made due to lack of data.Most of the works for plant tissue modelling have made use of sizers (see e.g.[24,25,48,88]), but there are several models that use timers as well [19,49].Although a timer will not allow cell size homeostasis over time and might trigger cell instabilities in the growing tissue [26], it can work as a mechanism for determining cell division times in a differentiating tissue with a countable number of cell cycles before cell division arrests.
In plant tissues, cell-to-cell interactions underpin the transport of key regulators and hormones that are fundamental at the tis-sue or at the whole plant scale.This includes, for instance, the directed transport of the phytohormone auxin.There are different theoretical proposals of how auxin is transported [89][90][91][92], and current studies are still evaluating through a combination of experiments and modelling whether such different competing existing models can explain patterning arising in different tissues [87,93,94].Other diffusible factors including cytokinins, proteins like the stem cell activator WUSCHEL, and microRNAs seem to be fundamental for patterning [95].Hence, in the study of certain patterning process, diffusible and directed transport between cells may need to be included into the models.
Modelling morphogenetic processes as a whole requires the simultaneous integration of growth, cell division, and signalling.An example where growth, cell division and intercellular signalling are important is the stem cell regulation in shoot meristems.In the SAM, a central pool of cells grows exponentially and the cells undergo division while the stem cell marker CLAVATA3 and WUSCHEL expression domains are continuously maintained.In this context, deterministic models have been able to predict the different expression domains in models including cell division and can also explain the variability seen following variability in tissue size [96].Ideally, a good (i.e.useful) model should at a first instance be sufficient to robustly describe the already existing data, and then, it should have some predictive power [14] (Fig 4).To achieve the first stage, one should define equivalent descriptors or observables in both experimental and modelled datasets that can easily be compared.A possible descriptor could be the expression pattern of a fluorescent reporter, or the time series of a certain variable of interest (Fig. 4).A more quantitative comparison can be performed by comparing experimental and theoretical concentration histograms and statistical properties of the different datasets.Classification analyses applied to both experimental and theoretical datasets can also bring quantitative ways to execute these comparisons (Fig. 4A-B, cf.Fig. 2E).Ideally, models should present explorations in the parameter space to establish whether the model can recapitulate the experimental data in a significant region of the parameter space, either by investigating large parameter regions [49,71,97,98], or if the model is more complex, by performing multiple optimisation runs [74,96].This will determine whether the model hypotheses are sufficient to robustly describe the experimental data.Still, several competing models can explain specific data sets equally well, and further tests of the model are required.Hence, in a second and more difficult stage, the model should be able to predict the outcome of experimental perturbations.Examples from the models we have discussed include the ability to predict an adaptive scaling of the stem cell domain to the size of the meristem, later confirmed in experiments [96], and a weak feedback component in the ATML1 sepal model that could later be verified in experiments [49].Other examples are the predicted growth patterns combining cell size and molecular input in the epidermal root tip (Fig. 4C-D) [37], and the prediction of auxin levels given the distributions of auxin transporters extracted from experiments (Fig. 4E-F) [82].Note that sometimes the theoretical parameter explorations can provide predictions of key experimental perturbations.For instance, in a recent study, a mathematical model predicted that the modulation of auxin influx transport would drive a change of vascular spacing in the Arabidopsis shoot, and this was experimentally corroborated in influx mutants [93].
As described, several examples of models of tissue morphogenesis exist where a robust simulated behaviour can explain data, and where novel predictions have come from the models.Still, the use of single parameter value explorations are dominating, in part due to the complexity of running models combining molecular regulation with growth.Hence, more efforts will have to be made in this direction.Also, comparisons between data and experiments can be further developed, and probably can be improved with the use of machine learning techniques.

Conclusion and outlook
In recent years, several efforts have been made for understanding morphogenesis in plants in a more quantitative manner.Many advances have been produced in both experimental and theoretical methods and applications, and these advances are enabling us to see the potential of cell-centred approaches to better understand tissue morphogenesis.Still, we need to further enhance and promote multidisciplinary efforts to get tangible and significant advances in the field of Computational Morphodynamics in plants.It will be essential to improve the ability to build on and directly compare published experiments and models from several groups to generate a coherent quantitative understanding of plant development.
Alternative experimental approaches, such as the use of plant cell strains, might still have an additional value to shed light into different puzzling phenomena that cannot be disentangled in planta.Indeed, single cell approaches, through the use of plant cell strains, have already brought very valuable knowledge into fundamental questions, including how auxin is transported in plants, and in how cell shape and cell polarity are generated [99].
A great challenge is to further automate experimental pipelines that can offer quantitative information in a systematic and highthroughput manner.The development of open source and userfriendly platforms for quantitative image analysis, and the creation of standards between such platforms for elaborating more complex pipelines will be pivotal for the establishment and consolidation of the Computational Morphodynamics field.Furthermore, the development of models with predictive power, close to the experimental data, will be key for a deeper understanding of plant morphogenesis.Such models might need to incorporate and integrate several mechanisms that can be present in different development scenarios, such as integration of mechanics and molecular regulation, dynamic stochasticity, alternative cell division rules, or even time-dependent parameters [66,[100][101][102], to emulate the evolving and the adaptable nature of regulatory networks in plants.
Finally, the success of the Computational Morphodynamics field will rely on keeping on bridging the gaps between experimental and theoretical approaches (Figs. 1 and 4) through communication between experimentalists and theorists, and through the development of new interdisciplinary scientists and laboratories.An integrative approach iterating between image collection, signal quantification, quantitative analysis and computational modelling.For each step, representative examples for methods and techniques are given and further detailed in the text.Note, this is not meant to be an exhaustive list of all available possibilities.

Figure captions
Figure 2.-Quantitative image analysis at the single cell level.(A) Image analysis pipeline for 3D nuclear segmentation from confocal images of developing sepals in a plant expressing the mCitrine-ATML1 fluorescent reporter [49].(B) 3D visualization of a shoot apical meristem (SAM) segmented using the watershed algorithm as implemented by Fernandez et al (2010) [45].(C) Single cell quantification of gene expression in the Arabidopsis root [82].Left: DII-VENUS (yellow) and cell geometries given by propidium iodide staining (red).Right: measured DII-VENUS levels quantified from the raw intensity image.(D) Single cell tracking of a SAM, using the block matching algorithm [18].Red lines represent newly formed walls within a period of 24 hours of meristematic growth.(E) Classification analysis of live imaging mCitrine-ATML1 expression data from sepal growth experiments (Fig. 2A) [49].Top left: ATML1 concentration time course of a small cell lineage.Top right: ATML1 concentration time course of a giant cell lineage.In both cases, coloured circles represent ploidy: yellow -2C; blue -4C; red -8C and above, and the grey lines show all cell linages.Bottom left: each circle represents the ATML1 concentration maximum recorded in 4C cells (i.e. during the G2 stage of the cell cycle) separately for small and giant cell lineages.Bottom right: performance of classification of cells (as either small or giant) based exclusively on ATML1 concentration maxima in 4C cells evaluated by area under the ROC curve (AUC).The red line is the ROC curve; diagonal dashed line represents AUC = 0.5 for comparison.An AUC of 0.8 suggests a good classifier (AUC = 0.5 corresponds to a random classifier; AUC =1 to perfect classification) and led to the hypothesis that a threshold-based mechanism for cell fate decision is at play involving a combination of high ATML1 concentration during a specific stage of the cell cycle.Horizontal black dashed lines in the top panels as well as bottom left panel represent the inferred ideal ATML1 concentration threshold.(F) quantification of topological features of Arabidopsis hypocotyls in different ecotypes [69].Top panel: tissue segmentation meshes with single cell segmentations; heatmap represents scale of betweenness centrality values.Bottom panel: virtual cross and longitudinal sections of extracted cellular networks for each ecotype.Heatmap represents scale of edge betweenness centrality values.Panels A,E have been extracted from [49], panel B from [45], panel C from [82], panel D from [18] and panel F from [69] 1 .Figure 3.-Modelling cell fate decisions in a growing tissue.(Left) Cartoon illustrating a growing tissue.A regulatory network is represented just in the central cell i for simplicity, in which a gene x autoactivates itself and activates a downstream target y.(Right) Possible schematic modelling pipeline.Top: Regulatory networks can be modelled either deterministically or stochastically, depending on the nature of the studied problem.In this scheme, a chemical Langevin equation is shown describing the x and y gene dynamics in cell i.Vi(t) is the volume of the cell i, f(xi) and g(xi) are generic functions describing the autoregulation of gene xi and the activation of yi by xi , respectively.ξz,i is a Gaussian random number with zero mean fulfilling 〈ξzi(t) ξz'j(t') 〉=δ(tt')δzz'δij , where i and j are cell indices, z and z' the modelled variables (x and y), δzz' and δij are Kronecker deltas and δ(t-t') is the Dirac delta [49].Tissue growth can be implemented by displacing the vertices out the tissue centre of mass, e.g., by following the proposed growth rule, where rj,k refers to the k-th coordinate of the j-th vertex in the tissue, and αk is the exponential growing rate along k-th direction.On top of the exponential growth, other mechanical cues can be implemented.In this case, we exemplify the mechanical contribution using a strain energy, S, which is a function of the strain (ε) and stress (σ) tensors.The integration of the growth rule explained above with the mechanical cues results into a more complex tissue growth, ui(t), where i refers to the i-th cell.This resulting growth needs to be taken into account when computing the rate of change of each concentration variable over time.When the resulting growth is more complex than an exponential, its dilution effects can be easily computed given the volumetric changes after each integration step.Finally, different cell division rules can be taken into account to determine the moment at which a cell divides and how the cell division plane will be positioned.Willis et al (2016) showed that the shoot apical meristem (SAM) divides using a rule that is in between the cell adder and sizer paradigms.The division plane can be set with the use of different algorithms [24].

Figure 1 .
Figure 1.-Typical Computational Morphodynamics workflow.An integrative approach iterating between image collection, signal quantification, quantitative analysis and computational modelling.For each step, representative examples for methods and techniques are given and further detailed in the text.Note, this is not meant to be an exhaustive list of all available possibilities.

Figure 4 .
Figure 4.-Bridging simulated data with experimental data.(A-B) Giant cell formation as a case study[49].(A) (Left) Scanning electron microscope image of a wild-type adult sepal, where giant cells are coloured in red.Scale bar: 100 µm.(Right) Simulation example of a developing sepal emulating the experimental data (cf.Fig.2E).The tissue grows anisotropically, while the ATML1 transcription factor dynamically fluctuates.Cells having ATML1 above a certain threshold when at the 4C stage (i.e.being in G2) are likely to endoreduplicate and become a giant cell.AU refers to arbitraty units.(B) Data analysis of simulation shown in (A) (cf.Fig.2E).Data analysis is performed on lower time resolution than the simulated time resolution in order to facilitate comparisons with the experimental data.(Top) Example of ATML1 simulated time-courses for a (top left) normally dividing cell and (top right) a cell becoming a giant cell.Colour refers to the ploidy of the cells as in panel (A) on the right.The red dashed horizontal lines refer to the soft threshold for giant cell fate commitment.(Bottom left) Spread plot showing the ATML1 maxima at 4C and the predicted ATML1 threshold.(Bottom right) ROC curve.(C-D) Mechanical study of the Arabidopsis embryo during seed germination [37].(C) Experimental data showing the relative cell expansion in the seed.(D) Mechanical growth simulations through finite element methods.Cell colours in Top panels in C-D show relative cell expansions.Growth rates depend on a combination of gene expression input and cell sizes [37].(E-F) Auxin patterning in the Arabidopsis root [82].(E) Auxin transporter localisations extracted from data.(Left) PIN auxin efflux exporters are shown in green.(Right) AUX1 (green and purple), LAX2 (blue), and LAX3 (purple) auxin influx importers.(Right) Scale bars = 50 µm.(F) Simulated DII-Venus levels (compare it with experimental data shown in Fig 2C) given the auxin transporter localisations shown in (E).