Virtual Plant Tissue: Building Blocks for Next-Generation Plant Growth Simulation

Motivation: Computational modeling of plant developmental processes is becoming increasingly important. Cellular resolution plant tissue simulators have been developed, yet they are typically describing physiological processes in an isolated way, strongly delimited in space and time. Results: With plant systems biology moving toward an integrative perspective on development we have built the Virtual Plant Tissue (VPTissue) package to couple functional modules or models in the same framework and across different frameworks. Multiple levels of model integration and coordination enable combining existing and new models from different sources, with diverse options in terms of input/output. Besides the core simulator the toolset also comprises a tissue editor for manipulating tissue geometry and cell, wall, and node attributes in an interactive manner. A parameter exploration tool is available to study parameter dependence of simulation results by distributing calculations over multiple systems. Availability: Virtual Plant Tissue is available as open source (EUPL license) on Bitbucket (https://bitbucket.org/vptissue/vptissue). The project has a website https://vptissue.bitbucket.io.


1.
VPTissue -VPTissue coupling: models  All model specific code (components and input files) can be found in src/main/models/Default/ (the files referring to TestCoupling, TestCoupling_I, and TestCoupling_II) TestCoupling_I and TestCoupling_II are static models so there are no Monte Carlo equilibration and growth or division rules required. Arbitrary units of time (T) and length (L) were used for the relevant processes. The simulation time step is 60 time units.
Intracellular reactions are described as: The left hand side of the equation represents the change in the concentration (number of molecules per area unit, N/L 2 ) of chemical 0 per time unit. A is the cell size (here area, L 2 ). The production rate is expressed on a per cell basis (N/T). Only the TestCoupling_I model produces chemical 0 ( 0 , prod k = 100 N/T).
To describe cell to cell transport we use the expression J from cell 2 to cell 1, with P a permeability constant, i c the chemical concentration of cell i. Intracellular diffusion is neglected in the models.
These equations describe the change in the chemical concentrations of each cell based on the summation of all passive transport flows over all its walls. These flow terms depend on a permeability constant of the chemical i i P (arbitrary high values were given that lead to quasi uniform concentrations of chemicals over the tissues: 0 P = 300, 30 L/T and 1 P =1000, 100 L/T for TestCoupling_I and TestCoupling_II, resp.), the concentration difference over the wall (   i chem  ), multiplied with the wall segment length l . The summed flow terms are divided by the cell's area to determine the chemical concentration change. To demonstrate flexibility of the coupling algorithm the identity of the transported chemical was defined differently in TestCoupling_I (chemical 0), and TestCoupling_II (chemical 1).
The coupling is defined in the input file of the TestCoupling model which identifies the two coupled models, the type of coupler (ExchangeCoupler), an arbitrary direction for the coupling (from TestCoupling_I to TestCoupling_II) the identity of the transported chemical in both models and the coupled cells (cell 28 to cell 0, cell 29 to cell 1, cell 30 to cell 2, cell 31 to cell 3). Each of the cell to cell couplings also has an associated transport constant ('diffusion'). The same equations are used as for transport in the individual tissues, but now using coupled cell concentrations, which are fixed during the coupling time interval.

2.
VPTissuepython (PyPTS): models 'leaf ' and 'root' (Fig. 6) The following equations describe the chemicals' concentration (number of molecules per unit area, N/L 2 ) dynamics through production, degradation, and transport in the individual models as well as transport between the models. As above (1) the units are arbitrary and the same expression for passive transport is used. Unlike for the previous models, the chemical production for these models is not on a per cell basis, but on a per area basis (this enables a smooth increase of chemical production with root growth). The time step is 1 time unit (T).
Kinetic constants for the leaf model:

 
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 VPTissue type of time evolution scheme was used (as in Figure 1), with node movement restricted to the y axis (mc_move_generator: 'directed_uniform').
Cells divide if their size is above 400 L 2 (with +/-10 % uniform noise added to avoid synchronous cell division).
The model specific code (components and input files) with a more complete listing of all parameters can be found in src/main/swig_sim/py_WrapperModel for the leaf model and in src/main/models/Default/ (files referring to WrapperModel). The boundary condition coupling is implemented in a different way than above with the lower cell row of the leaf model ( The colouring represents increasing morphogen concentrations by increasing saturation of the yellow versus the red colour. Uncoloured cells have a morphogen concentration below an expansion threshold. The snapshot in (A) corresponds to 668 simulation steps. Cells of the root model contain an arrow representing the PIN transporter polarization. Model simulations with over 1000 cells like the leaf simulation here (snapshot with 1058 cells after 67 simulation steps) can take less than 10 minutes with Virtual Plant Tissue (depending on model and on a standard desktop computer with an Intel Core i5-5300U CPU @ 2.30GHz, with 8 GB RAM). The models are part of the Virtual Plant Tissue distribution (for details see user manual). More precise run time comparisons for specific models and cpu's can be found in Suppl. Table 1.

Supplementary Figure 3.
Figures A-D illustrate the fundamental differences between the Hamiltonians, in a very simple and direct way, based on the mechanical equilibration of the Geometric model with one small cell and one larger cell which are connected by a wall (A; start of simulation). Only the internal nodes in the connecting wall can move during simulation. After 150 time steps with the original Hamiltonian ('PlainGC') the larger cell already strongly penetrates the smaller cell (B). The same is true for the 'ElasticWall' Hamiltonian (C; the effect is not very pronounced yet due to stronger elastic counterforces). This is no longer the case for the 'ModifiedGC' (D) and the 'Maxwell' (E) Hamiltonians where the smaller cell converges to a stable convex appearance as seen experimentally [Corson et al., 2009]. The more subtle shape difference in D shows that, unlike for the ModifiedGC, through the more realistic biophysical properties of this Hamiltonian it is possible to fine-tune the convexity of the small cell.