Fundamentally, a plant tissue is described in Virtual Plant Tissue as a mesh consisting of polygonal cells that consist of nodes (vertices) connected by edges. The cells have attributes (geometric, chemical, and others) and neighboring cells share an apoplast segment or “wall” (consisting of one or more edges). The walls in fact represent two physiologically distinct cell wall (and plasma membrane) domains on either side of the middle lamella and therefore have neighbor cell specific attributes as well as others (Supplementary Figure 1). Thanks to the “transporter” array attribute, PIN-FORMED (PIN; Krecek et al., 2009) transporter levels can be assigned to the cell wall segments and therefore, like VirtualLeaf (Merks et al., 2007; van Mourik et al., 2012), Virtual Plant Tissue is suitable to study auxin transport and PIN polarization dynamics (the “AuxinGrowth” model is derived from Merks et al., 2007). However, unlike other modeling approaches such as cellular Potts models (Grieneisen et al., 2007) or vertex-based modeling with cell partitioning (Abley et al., 2013) intracellular hormone concentration gradients cannot currently be represented. On the other hand the latter types of models are typically associated with a high computational cost and have other issues (for instance cell sliding in lattice based models). Virtual Plant Tissue can efficiently simulate large tissue structures (Supplementary Table 1, Supplementary Figure 2), however, some implicit assumptions have to be considered. An implicit assumption is that intracellular gradients are small, which for auxin may not be the case for cells as long as 100 μm (Kramer, 2008). Another inherent assumption is that lateral diffusion in the apoplast is negligible which in principle needs to be assessed case by case.

The state of the tissue at a specific time is described in a single XML input file or “leaf file” (in accordance with Merks et al., 2011; the leaf files are, however, not cross-compatible). An input file contains model specific parameters (including names for the components or modules that define the model), but also lists all cells, nodes and walls together with their unique attributes (cf. user manual: Supplementary File 1). Fundamentally, the dynamic processes during a simulation time step are separated into a phase with slow biochemical and biophysical processes governed by ordinary differential equations and algebraic rules, and a phase with instantaneous (elastic) wall mechanical processes, governed by a Monte Carlo sampling strategy to approximate the equilibrium state (Figure 1). This type of time-scale “separation” and the use of a generalized energy function or Hamiltonian is an efficient and established modeling approach (Wabnik et al., 2010; De Rybel et al., 2014). Since this only excludes simulating very fast (elastic) mechanical changes (taking <1 min: see for instance Proseus et al., 1999), it does not strongly limit the scope of plant biological questions to be studied. The ordering and composition of the phases is customizable as explained below.

The Virtual Plant Tissue code structure has been designed for a modular model building and simulation approach at three principal levels: (i) combining functional units within a model, (ii) coupling models within Virtual Plant Tissue, and (iii) coupling with external models of other frameworks (Figure 1). In Virtual Plant Tissue biological as well as algorithmic concepts/modules have been molded into classes (Table 1, Supplementary File 1). Importantly, the regulation of the time evolution scheme, defining type and sequence of computational processes (cell division, wall yielding, hormone transport, …), is encapsulated in a separate time evolver class. Generally, a great flexibility to define model simulations is available to the user by editing parameters listed in the input XML file. Importantly, all parameters (including model/module selection and time evolution algorithm) can be changed dynamically, with every change taking effect at the start of the time step following that change. This can be done either interactively with the “gui” or as the result of an automated computation, i.e., when certain conditions are met such as surpassing a threshold value of number of cells or simulated time. One could for instance define a component of components, for example a cell division (cell_split) module which selects either symmetric or asymmetric division (both defined in separate modules) based on cell type (pavement cell or meristemoid cell).

The core simulator has been isolated into a single class which enables multiple instances of simulations to be run in parallel via an internal interface. This represents a higher level of modularity linking internal models/simulations. A Coupler type class defines how multiple simulations can communicate. An example coupler class is available that maximally preserves the respective models' individual character by restriction of information exchange through modified boundary conditions (ExchangeCoupler class (Figure 2). With a dedicated parameter in the VPTissue input file the coupling frequency can be optimized (Figure 5). This is a necessary feature since the coupling step (time slice) needs to be small enough to ensure mass conservation across the boundary. In the example of Figure 5A one of two basic coupled models only produces a chemical (model details in Supplementary Text 1.1). The accumulation of that chemical in the uncoupled case is only closely approximated if the coupling constant is high enough (leading to intensive updating of boundary condition). As expected increasing the permeability constant (“diffusion” in the xml file) for chemical exchange between models leads to nearly equal partitioning of the chemical (Figure 5B).

A SWIG adapter has been implemented to define a clean external interface opening up the possibility for coupling multiple models implemented on different platforms. The interface has methods that are also defined in a wrapper class that allows the Virtual Plant Tissue core simulator to be accessed from an external program that can be written in a different language (Python or Java) via native wrapper objects (Figure 4). As an example a model defined in Python (with the PyPTS package,) describing a leaf was coupled with a simple root model defined in Virtual Plant Tissue (Figure 6, Supplementary Video 1, model details in Supplementary Text 1.2). In this particular case the Python model is supplemented with instructions to coordinate the coupling and simulation events (a separate driver or organizer script would work equally well). The (geometrically static) Python leaf model accumulates chemical “0,” which is transported to the Virtual Plant Tissue root model. The root responds above a threshold by growing and increasing chemical “1” production, which will be transported to the leaf and eventually will shut down chemical “0” production and subsequently root growth.

Concrete biological questions through the described approaches for coupled simulations are diverse. For instance: how does inter-organ transport of assimilates and hormones affect leaf growth (e.g., Bhalerao et al., 2002; Ljung et al., 2015; Notaguchi and Okamoto, 2015)? This could be tackled by coupling an architectural plant model (functional-structural plant model; Vos et al., 2010) of the shoot with a vertex-based model of leaf growth in our framework (a starting leaf model is already available). The most likely candidates for exchange through the external (wrapper) interface are sugars and hormones with the fluxes depending on the developmental state of the respective organs. A direct internal coupling of VPTissue based root and shoot models is feasible too (with minerals and hormones as candidates for exchange). Another possibility for internal coupling is to study the influence of plant architecture (branching) on polar auxin transport (Bennett et al., 2016). Stem segments could be coupled in a cell-based way to exchange auxin. Some of these ideas are illustrated in Figure 7.

Despite the extensive efforts to make the framework readable and minimizing dependencies in the code, developers adopting the platform will still need to spend some time familiarizing with it. Basic tasks to master are adding models, model components and model attributes to the code base. Adding a model requires defining an input file as well as some preferences files (in src/main/models/'model family'/resources). One can copy an existing model as the template, use an xml editor or use the Virtual Plant Tissue Editor (see below) for that. New model component files need to be added to src/main/models/'model family'/components/'component type'. Each new component also needs to be added to the local factories.h and factories.cpp files and to the CMakeLists.txt file one level higher. A convenient feature made possible through the Virtual Plant Tissue code modularity is the possibility to define a (meta-)component that uses other model components from the same component type. For instance, one could define a leaf growth model that contains different cell types such as pavement or meristemoid cells. A cell division meta-component could then simply refer to separate dedicated pavement and meristemoid division modules without the need to combine them and duplicate a lot of code. The steps for adding an attribute depend on the type (belonging to simulation, cell, wall, tissue, or node). Detailed help on programming with Virtual Plant Tissue can be found in the user manual (Chapter 7).

Diverse models are provided in the current Virtual Plant Tissue distribution (some originally from VirtualLeaf). Several of them are meant for demonstration or as a starting point for building more advanced models, for instance to study phyllotactic patterning (model based on Smith et al., 2006) and leaf venation (with Meinhardt and AuxinGrowth models). Other models correspond exactly to what was published, including root and leaf models. They can also be used to study regulatory networks involved in leaf and root growth. Some model simulation snapshots can be seen in Supplementary Figure 2. The main components defining a model are listed in Table 1 and can be divided into biological and Monte Carlo mechanics modules. Thanks to the modularity of the Virtual Plant Tissue code a new component type can be easily added and used by different models or model versions. For the biological component types cell_chemistry, wall_chemistry and cell2cell_transport the time organization within each simulation time step is determined by a differential equation solver. On the other hand, components like cell_housekeep, cell_split and cell_daughters contain rules that are evaluated only once per simulation step. For cell division (cell_split) the options currently available are still the same as in VirtualLeaf: cells divide into equal daughter cell areas based on a pre-specified orientation of the division axis or a division axis perpendicular to the cell's axis of inertia (Merks et al., 2011). For the mechanical components different options are now available. The original Hamiltonian of VirtualLeaf is now called “PlainGC” and uses the difference of absolute target cell areas and actual cell areas (Merks et al., 2011). In the “ModifiedGC” the relative differences of the cell areas are used which avoids terms of larger cells to dominate the Hamiltonian (Supplementary Figure 3). Original:H=λA∑i(a(i)-AT(i))2+λM∑j(l(j)-LT(j))2,Modified:H=λA∑i(a(i)-AT(i)a(i))2+λM∑j(l(j)-LT(j))2, where indices i and j sum over all cells and polygon edges, respectively, λ_A is a parameter setting the cells' resistance to compression or expansion, and λ_M is a spring constant. A_T is the cell's target area, L_T the wall element target length. The target length is constant and the same for all wall elements (edges) and new nodes can be inserted as soon as an edge reaches a critical length (“target_node_distance”). In the Hamiltonian of the “ElasticWall” model the edge length constraint is replaced by an elastic energy term and each wall has its individual wall strength and its variable rest length (specified in the xml data file). This provides more flexibility to describe wall mechanics. The “Maxwell” Hamiltonian additionally replaces the target area constraint by a turgor pressure term in analogy with Corson et al. (2009): Maxwell:H=-∑cellsPS(i)+Eh2∑wallsl0(j)(l(j)l0(j)-1)2. With the elastic modulus, h the uniform wall thickness, l(j) and l₀(j) the actual and rest lengths of wall j, respectively, the turgor pressure, and S(i) the area of cell i.

Every time step each cell's solute quantity and each wall's rest length are updated through the corresponding cell_housekeep module. For the ElasticWall model irreversible wall extension happens by increasing wall rest length by a constant fractional amount if the actual wall length exceeds a threshold (for instance if 50% higher than the rest length). In the Maxwell model the wall rest length and solute concentration are updated according to the solution of differential equations (Equations 2, 8) described in Corson et al. (2009) as applied to one Virtual Plant Tissue simulation time step. The repertoire of mechanical modules is useful to investigate the influence of wall mechanics on organ growth. A more concrete question is how differences in wall mechanical properties and osmoregulation of different sections of the maize leaf affect growth under abiotic stress conditions such as drought (e.g., Dzhurakhalov et al., 2015b). The Maxwell module is particularly suitable since, instead of the artificial representation of wall yielding and turgor in VirtualLeaf, real biophysical properties are used such as elasticity and viscosity.

Whereas, Virtual Plant Tissue offers the speed and flexibility to develop more advanced models for cell wall mechanics, it currently does not allow describing intracellular interactions which potentially play a role in determining cell shape such as those involving the cytoskeleton (Sampathkumar et al., 2014). To describe such processes at high resolution Virtual Plant Tissue is less suitable than finite element methods (e.g., Yanagisawa et al., 2015) or a three-dimensional modeling framework (Boudon et al., 2015).

In line with the modular simulation set-up diverse new algorithmic options are available, like choices for the ODE solver (www.boost.org), random number generators and distributions, and Monte Carlo energy evaluation criteria. It is for instance important for a (stochastic) modeling framework to evaluate the influence of those choices on simulation output and to ensure convergence of the Monte Carlo equilibration (Dzhurakhalov et al., 2015a). For an overview of algorithmic options the reader is referred to the user manual (Chapter 3). In the end all selected model components are organized in a time evolution scheme (time_evolver class). Virtual Plant Tissue provides different readymade choices to the user. For instance selecting the “VLeaf” evolver instead of the “VPTissue” evolver results in each simulation time step terminating with the reaction and transport steps instead of the fast elastic equilibration step.

Simulations are organized into workspaces which consist of projects (directories) comprising the initial data file, simulation output preferences, and accumulated output files. Besides running Virtual Plant Tissue via the command line a graphical user interface is available for users. Figure 8 shows a Virtual Plant Tissue screen shot open at a workspace with several projects (left panel). The top right panel is used to access the workspace preferences determining features such as which I/O viewers need to be enabled, which color scheme should be used, etc. The Parameters panel allows viewing and editing all configuration parameters from the simulation. The Project Preferences panel enables to overrule workspace preferences for a particular project. The bottom panel provides a running log of a project that is open. By virtue of a rigorous Model-View-Controller design it is possible to attach and remove multiple viewers during a running simulation to economize on computational resources. Figure 9 shows a screenshot of the Qt viewer for the Geometric project in Figure 8.

The simulator itself has functionality for converting between XML, compressed XML or HDF5 data formats (Supplementary File 1) and for post-processing simulation output in various graphic (PLY, pdf, png, …) or text (csv) formats. The HDF5 format allows data arrays (practically unlimited in size) to be easily accessed, exchanged and analyzed (a plugin for the Paraview visualization software is included in the source code, a dedicated VPTissue HDF5 file reader for Python, is available on https://pypi.python.org/pypi/PyPTS/0.2.4).

The Virtual Plant Tissue package comes with a toolset comprising a graphical editor (simPT_editor) for XML files. The Tissue Editor is a graphical editor for the Virtual Plant Tissue mesh geometry, the cell, wall, and node attributes and the model and simulation parameters (Figure 10). The application constructs, reads and writes a full XML file, including simulation parameters and mesh data. It is possible to load an image that can be used as a template to draw cell meshed based on microscopic images of plant tissue. Detailed information on the graphical interface can be found in the user manual (src/doc/latex_user_man/UserManual.pdf, Chapter 3).

Virtual Plant Tissue also comes with a parameter exploration tool (simPT_parex) which allows running and monitoring a parameter sweep calculation (with different sampling strategies) on a compute server. After connecting to a server one can start a new exploration or get details on a running task. Several types of exploration can be started: in sweep based explorations one parameter is varied based on a range of values (e.g., Supplementary Figure 4), template based exploration allows varying several parameters at once with a csv file that specifies the parameter combinations. More detailed information can be found in the user manual (src/doc/latex_user_man/UserManual.pdf, Chapter 4, Data Sheet 2).

The build and installation process is tailored to rapid identification of compilation/installation/run-time errors with a high level of abstraction and platform independence allowing easy compilation on Linux, Windows, and MacOS systems. Code documentation is partly automated [the Application Programming Interface (API) documentation in the reference manual is generated via doxygen (www.doxygen.org)].

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.