OpenSimRoot: widening the scope and application of root architectural models

Summary OpenSimRoot is an open‐source, functional–structural plant model and mathematical description of root growth and function. We describe OpenSimRoot and its functionality to broaden the benefits of root modeling to the plant science community. OpenSimRoot is an extended version of simroot, established to simulate root system architecture, nutrient acquisition and plant growth. OpenSimRoot has a plugin, modular infrastructure, coupling single plant and crop stands to soil nutrient and water transport models. It estimates the value of root traits for water and nutrient acquisition in environments and plant species. The flexible OpenSimRoot design allows upscaling from root anatomy to plant community to estimate the following: resource costs of developmental and anatomical traits; trait synergisms; and (interspecies) root competition. OpenSimRoot can model three‐dimensional images from magnetic resonance imaging (MRI) and X‐ray computed tomography (CT) of roots in soil. New modules include: soil water‐dependent water uptake and xylem flow; tiller formation; evapotranspiration; simultaneous simulation of mobile solutes; mesh refinement; and root growth plasticity. OpenSimRoot integrates plant phenotypic data with environmental metadata to support experimental designs and to gain a mechanistic understanding at system scales.


Note S2: Command line interface (CLI) of OpenSimRoot: How to run and use the model
OpenSimRoot has a command line interface, which means that you operate the model from a terminal using commands, not with a graphical interface and the mouse.
Step 1: Open a terminal (under windows 10 you may use the program named CMD) Step 2: Go to the folder where you want to run the model, use the command cd to navigate, for example: cd MyRunFolder Step 3: We assume that the folder contains the OpenSimRoot executable. With the "ls" command you can list all folders (or on windows the command is "dir"). Here we see that my folder contains the executable OpenSimRoot (conveniently made green, as it is executable) and a XML input file.
Step 4: OpenSimRoot has a small build in help which we we can run by typing ./OpenSimRoot -h (on windows you do not type the path "./" in front of the executable).
The help shows how to run OpenSimRoot, and gives you some options and their explanation.
Step 5: Like the help shows, running the model is done by appending the input file: ./OpenSimRoot SimpleCropModel.xml Again with ls (dir) you can list the filer, the model created two new files, one containing warnings, one containing the simulation results.
Step 6: The results of the simulation are in the tabled_output.tab file which can be viewed with any program that opens text files. Here we simply show the first lines with the command head: The file contains a header in the first row, and 6 columns listing the name of the state variable, the time, the value, the rate of change of that state variable (if simulated), the unit of the state variable, and the path in the hierarchy to this state variable.
Real time hours: minutes : seconds Results file One warning Simulation time

Command
Step 7: The tabled_output.tab file is also easily imported into a spreadsheet program. By enabling the auto filter and selecting leafArea, we can easily create a plot. 8. The same can be achieved in R using this script: d<-read.table("tabled_output.tab",header=T) f=d$name=="leafArea" plot (value~time,data=d[f,],ylab=~"leaf area (cm"^2*")", xlab="time (d.a.g)") step 9: Editing the input file can be done with any text file editor. Here I opened the file with the command nano tabled_output.tab and the result is an xml formatted file in which we can change the numbers, save and rerun. In white you see the numbers, and scrolling to the bottom you would see more.
step 10: You see several functions listed that are used to simulate a state variable. To get a list of all functions that are included in your OpenSimRoot version use the command OpenSimRoot -l. This will list all plugins that are included with OpenSimRoot.

Note S3: Class hierarchy of OpenSimRoot code
This document lists the class hierarchy for the most important classes in OpenSimRoot.

Minimodels
Minimodels are of type SimulaBase and encapsulate one time and location dependent state variable. The inhertance diagram for all SimulaX classes is given in Figure S3.1. • SimulaExternal provides a mechanism for encapsulating other dynamic simulation models. • SimulaPoint simulates a point and its movement through space.
• SimulaVariable simulates a value and change over time using numerical integration.
• SimulaGrid simulates a static, 3D field using a list of Coordinates with values and a 3D interpolation algorithm • SimulaLink simply bridges to another minimodel in the hierarchy of minimodels.
• SimulaStochastic draws numbers from a random number generator.

Inherited from DerivativeBase
Below is a list of all the plugins that directly, or indirectly, inherit from DerivativeBase and can be used by SimulaVariable, SimulaPoint or SimulaDerivative for computation.

List of plugins for simulating various processes
Note that these are a list of classes, as they appear in the code. Registration of the plugins may occur under different names. Inputfiles use the registered names, not the class names. Use OpenSimRoot -l to get that list. See also operation manual in Note S2.

Integration functions
The SimulaVariable and SimulaPoint classes use helper functions for integrating the result. Several integration methods have been implemented ( Figure S3.2). New integration functions can be added and registered, using the plugin framework, similar to the classes that inherit from DerivativeBase.

Object generators
Object generators are plugins that can be associated with any SimulaX object and update the list of children when a child is requested.

Note S4: Plugin example code
Here we give example code for a simple plugin and the code needed to register this plugin with OpenSimRoot. Once the code has been put into a text file, it can be compiled and linked to OpenSimRoot.

1) For new algorithms
//Class declaration. Class should inherit from DerivativeBase, have a constructor, and implements two virtual methods, getName() and calculate(). The example class presented here has two SimulaBase pointers as private members, which will be used to connect to the minimodels that simulate length and diameter of a root segment and to retrieve their values. //the constructor of our class. pSD is the pointer to the minimodel that uses the plugin for computation RootSegmentSurfaceArea::RootSegmentSurfaceArea(SimulaDynamic* pSD):DerivativeBase(pSD) { //We check that the user set the unit right pSD->checkUnit("cm2"); //We retrieve the pointers length=pSD->getSibling("rootSegmentLength","cm"); diameter=pSD->getSibling("rootDiameter","cm"); } //the computation void RootSegmentSurfaceArea::calculate(const Time &t,double &area){ //first we retrieve data double d,l; diameter->get(t,d); length->get(t,l); //second we compute area=l*d*PI; } //the name under which the plugin will be registered, make sure it is unique, use OpenSimRoot -l to see what names are already taken std::string RootSegmentSurfaceArea::getName()const{ return "rootSegmentSurfaceArea.v3"; } //Now we create a function for instantiating our class DerivativeBase * newInstantiationRootSegmentSurfaceArea(SimulaDynamic* const pSD){ return new RootSegmentSurfaceArea(pSD); } //And we register this plugin using a static instantiation of a class which guarantees that the constructor is when OpenSimRoot is started static class AutoRegisterMyNewPlugin { public: AutoRegisterMyNewPlugin() { //this line does the registration. Make sure you register under the same name as the getName() method returns. This important for the model dump being loadable again.

Note S5: Detailed description of the water and nutrient submodules
Watermodule Plant transpiration is simulated by OpenSimRoot, assuming that water availability is not limiting and stomatal conductance is constant. Transpiration and evaporation need to be separated within OpenSimRoot. Transpiration can be estimated from a fixed water use efficiency parameter (which simply links carbon fixation linearly to transpiration), or from the Penman-Monteith model, which computes evapotranspiration based on weather conditions (Penman, 1948;Monteith, 1964). When transpiration is calculated based on a water use efficiency parameter, the user needs to provide evaporation values; when the Penman-Monteith model is used, transpiration and evaporation are separated by OpenSimRoot solving the Penman-Monteith model twice, once for full crop cover, and once for a bare soil. Based on the percent light capture by the crop OpenSimRoot scales evaporation and the transpiration terms assuming evaporation is negligible and small under full crop cover (Leaf Area Index ~3). To simulate the soil hydrology, OpenSimRoot has a submodule that solves the Richards equation in three dimensions using finite element method (FEM) on a Cartesian grid. The soil water submodule is a simplified and modified C++ rewrite of the SWMS3D model, which is the basis of Hydrus and R-SWMS (Šimunek et al., 1995;Diamantopoulos et al., 2013). Certain exceptional circumstances such as drainage or water ponding at top soil, are excluded. The top boundary condition is a water flux that is the difference between precipitation and evaporation. Evaporation, as computed by the Penman-Monteith equation, is assumed to be potential evaporation (i.e. appropriate for wet soils), and assumed to be equal across the soil surface, shoot geometry is not simulated. Potential evaporation is scaled back to an actual evaporation by including a smooth scaling function which causes evaporation to decrease smoothly from potential, when the top soil is wet, to equal the soil conductivity when the soil is not able to sustain higher evaporation rates. If the top soil is not necessarily uniformly wet, actual evaporation will be non-uniform across the soil surface in OpenSimRoot. The water retention curve and soil hydraulic conductivity are computed using the van Genuchten and Mualem equations. The Richards equation can include a sink term, which in OpenSimRoot represents water uptake by roots (as described evaporation sink is handled as dynamic boundary condition). To do so we need to know 1) how much water is taken up by each root segment at a given moment in time, and 2) how that uptake is coupled to the FEM nodes of the grid on which the Richards equation is solved. Assuming that root uptake equals transpiration, i.e. we ignore temporal water storage in the plant, OpenSimRoot can either divide the water uptake of the whole root system by assuming each root segment contributes equally to uptake relative to its length (as in Hopmans, (Hopmans & Bristow, 2002)) or by solving the hyrdraulic architecture represented by a network model and using a circuit analogy likewise motivated by finite element theory (Alm et al., 1992;Doussan et al., 1998). The network model is novel in OpenSimRoot implemented to work with a growing root and used in the study of Schneider (Unpublished). This model requires axial and radial hydraulic conductivities for each root segment, which can be defined in the input files as a function of root age and class, and are scaled (i.e. normalized) with the inverse of the root segment length (axial), or the root segment surface area (radial). The coupling of the root model to the FEM model enables each root segment to have a soil water content at the root surface. The next step is to make sure that water uptake by the root system equals the transpiration which is achieved by changing the water potential at the root collar (top of the hypocotyl). Getting the root collar potential is a parabolic optimization function which is solved with a newton solver, typically in three steps. The water potential at the top of the hypocotyl is not allowed to drop below a given threshold. If the threshold is reached, OpenSimRoot assumes that water uptake is less than potential transpiration and will write a warning. Further simulation results might not be correct as currently no effects of drought on photosynthesis, leaf expansion etc have been implemented. However, the model should correctly deal with compensatory uptake of water when soil water distribution is heterogeneous. And this model can show water loss of roots while the same conductivity from xylem to soil is assumed. Mapping the root model to the FEM model is done based on a neighborhood search. All FEM nodes surrounding the root segment are considered. Sink terms, and local environment are computed based on inverse distance weighted average of the FEM nodes surrounding the root node. An alternative mapping algorithm, by which every FEM node is assigned with every root node has been implemented, in order to ignore root architecture completely in the water and nutrient uptake simulations. This was for example used in Postma and Lynch (2012) where it was concluded that the positioning of the root, that is root architecture, is necessary for simulating niche differentiation for nitrate uptake among maize, bean and squash plants, whereas if roots would be able to take up nutrients from everywhere in the soil, there would be no niche differentiation.

Nutrient module
OpenSimRoot has a nutrient module to simulate the uptake solutes, and in the new version theoretically simultaneously for various nutrients. This module was implemented to simulate the function of root architectural traits for nutrient uptake, and test tradeoffs for acquisition of different nutrients. Time dependent optimal and minimal nutrient content (µmol/g) have to be defined for leaves, stems and all root classes, for to be simulated solutes. These amounts are used to compute nutrient requirements of the plant, and compared to total uptake amounts, including initial seed reserves (for uptake see below). When uptake is less than demand, plant stress is assumed, with maximum stress being defined as uptake equal to minimal nutrient content (stress(uptake) = max( 0, min((uptake-minimal)/(optimal-minimal),1) ). Stress modifying impact functions can be defined for components such as leaf expansion rate, photosynthesis rates, respiration rates, and root elongation rates or secondary growth. Typically, they should be defined such that, when stress=0, growth ceases altogether. For example, by making the initial response of the shoot stronger than that of the root, the plant will decrease shoot to root ratios when nutrient deficient. Thus OpenSimRoot will move towards a functional equilibrium, although due to the inherent slow nature of growth, and the relative fast dynamics of other processes, this functional equilibrium might not be reached, and oscillatory behavior might occur (Postma & Lynch, 2011;Postma et al., 2014b). The current implementation assumes that internally, reallocation of nutrients is fast and perfect, such that all organs experience equal stress. This might be true for a nutrient like nitrogen, which typically causes chlorosis everywhere in the shoot, but might not be correct for other nutrients. The importance of simulation of nutrient redistribution in the plant still needs study, and would require implementation of a shoot architectural model in which the age and position of individual leafs is tracked. Nutrient uptake from soil to root is simulated independently of utilization of nutrients within the plant. Two options for simulation are provided: 1) The Barber-Cushman model and 2) a 3D FEM model. One is a C++ implementation of the original Barber-Cushman model with root hairs. The model is described as radial 1D PDE (Partial Differential Equation) which corresponds to the rhizosphere around the root. It assumes nutrient uptake to be described by a Michealis-Menten term, and the nutrient transport in the soil to be driven by convection (water flow) and diffusion. A buffer constant replaces a reaction term. The Barber-Cushman model is suitable for immobile nutrients like phosphorus. Phosphorus uptake causes steep gradients in concentrations around the root. These depletion zones are typically only 2-4 mm in diameter, and thereby would require a computationally unacceptably high resolution of the 3D finite element model (~0.1 mm resolution of a 1 m 3 soil pedon would result in 1e 12 elements or 8 petabytes to hold a single double precision array). Competition between roots is computed based on a local average root density which determines the outer boundary of the Barber-Cushman model. OpenSimRoot updates this boundary when new roots grow in the vicinity of other roots and corrects the initial nutrient concentration for new roots with the uptake of nutrients of older roots. Nevertheless, this handling of root competition is only acceptable when the overlap of depletion zones, which can be computed based on raster images of the root system, is relatively small. For crops, overlap in phosphorus depletion zones is typically below 20% because of its low mobility. Inter and intra root competition plays a much more important role in the uptake of mobile nutrients such as nitrate. Nitrate might form diffuse or no depletion zones around the root and for this reason is better simulated using a 3D FEM. SimRoot solved the convection-dispersion equation on the same FEM grid as the water transport is solved which can be restricting, OpenSimRoot alternatively can solve it on a refined grid, where the refinement factor is yet fixed to 2nd, 4th, 8th or the 16th of a reference grid. For each solute a new FEM model is instantiated and linked to the water model. The 3D FEM model for solute transport is coupled to the root systems using the same method as used for the hydraulic model, where the uptake of solutes by the root segments is based on Michaelis Menten kinetics, as in the Barber-Cushman model. Buffering and diffusion coefficients are dependent on the soil water content, and might thereby deviate from the constant coefficients used in the Barber-Cushman model. The effects must be considered when comparing the output of both models (Postma and Lynch, 2011). When simulating more than one solute, solutes do not influence each other directly in OpenSimRoot. Indirect effects occur through the influence of nutrient uptake on root growth. Each solute has a stress function to determine how each impacts, for example, photosynthesis. A user specified aggregation function determines the aggregate impact (Dathe et al., 2013). For example, Postma et al., (2014a) showed how the optimal lateral branching density in maize depends on the relative availability of phosphorus and nitrogen.

Note S6: Example of a simple OpenSimRoot input file
The XML below is an example of an OpenSimRoot input file that constructs a simple crop model, without any roots. All the SimulaX tags will instantiate a minimodel of the corresponding type, for example a constant (time independent parameter) is declared as <SimulaConstant ...>. Metadata for the minimodels, such as name and unit, are given in the attributes lists.
General rules for XML documents 1) The document has tags which are between brackets like <> 2) Tags correspond to minimodels in OpenSimRoot and therefore carry different names, such as SimulaBase, SimulaConstant, etc.
3) Tags need to be closed either by putting a / before the closing bracket, or if data is nested inside the tag with a corresponding closing tag which is recognized by </. For example <SimulaConstant></SimulaConstant> 4) Between opening and closing tags you will find data, and or declarations of minimodels which are at the next level in the hierarchy 5) Tags carry attributes which describe metadata. Attributes are always listed as attribute="something". In OpenSimRoot all tags have at least a name attribute. 6) An XML document is plain text and recognized by a special declaration at the top of the document. <?xml version="1.0" encoding="UTF-8"?> 7) XML documents can have stylesheets associated with them so the the browser knows how to render the document. Here we have <?xml-stylesheet type="text/xsl" href="tree-view.xsl"?> 8) Comments are between <!--and -->. 9) All XML documents of a document type tag. For OpenSimRoot the document type is declared as <SimulationModel></SimulationModel>. All other tags must be in between these tags.
Here follows an example input file. The comments in black give more explanation as to how a simple crop model is being constructed by this input file. Input files for full root architectural models can be found in the software repository on gitlab: https://gitlab.com/rootmodels/OpenSimRoot <?xml version="1.0" encoding="UTF-8"?> <?xml-stylesheet type="text/xsl" href="tree-view.xsl"?> <!--Copyright © 2016 Forschungszentrum Jülich GmbH All rights reserved.
Redistribution and use in source and binary forms, with or without modification, are permitted under the GNU General Public License v3 and provided that the following conditions are met: 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution. 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission. You should have received the GNU GENERAL PUBLIC LICENSE v3 with this file in license.txt but it can also be found at http://www.gnu.org/licenses/gpl-3.0.en.html --> <!--This XML constructs a simple, radiation use efficiency based crop model.
Roots and stems are only presented as Carbon (dry weight) pools Leaf dry weight is converted to leaf area based on specific leaf area (SLA) Leaf area is converted to light interception using an extinction coefficient. Light interception is converted to photosynthesis using radiation use efficiency (RUE). Photosynthesis is converted to structural carbon using a conversion factor (multiplier) which represents relative losses due to respiration Fixed allocation causes structural carbon to be divided over root, stem and leafs.
Behavior, is simple exponential growth for which RGR = SLA * C2Leafs * photosynthesis * multiplier However, as the light interception with increasing leaf area reaches an asymptote, the model will move towards linear growth.--> <SimulationModel> <!--SimulaBase is a simple container, that holds other SimulaX objects. SimulaBase is thus a minimodel that does not hold or simulate data. It should, like all mini models, have a name. So here we declare a container in which we are going to put all our plants. Inside it we put a container for our plant, named arbitrarily "myPlant".
--> <SimulaBase name="plants"> <SimulaBase name="myplant"> <!--Here follow three SimulaConstant declarations. SimulaConstant is a minimodel that holds time and space independent data of different types. Possible types are double, int, string, Coordinate. Besides the name attribute they must have a unit, and if the data is not a double, a type declaration.
A plant should be of a given species/genotype. The model will look for a parameter set in roottypeParameters with the corresponding type. Here we declare that we want to simulate a plant of type mySpecies --> <SimulaConstant name="plantType" type="string"> mySpecies </SimulaConstant> <!--The time that the plant is planted. 0. is at the start of the simulation. --> <SimulaConstant name="plantingTime" unit="day"> 0. </SimulaConstant> <!--Location in space where the seed is planted. --> <SimulaConstant name="plantPosition" type="Coordinate"> 0 -2 0 <!--Container that hold all the minimodels that will simulate shoot related parameters. The shoot and root are inside plantPosition, as OpenSimRoot works with a relative Coordinate system. We achieve that all coordinates that belong to our plant are relative to plantPosition. --> <SimulaBase name="shoot"> <!--Licht interception is simulated by the light interception module. SimulaDerivative declares a minimodel that wil use the lightInterception plugin to compute light interception. Attributes are name of what is being computed (name="lightInterception"), the unit of what is being computed (unit="umol/cm2/day"), and the plugin that should be used to compute it (function="lightInterception"). The plugin lightInterception requires leafAreaIndex and from the parameter section and extinctionCoefficient (kdf). Further it needs irradiation levels from the environmental section. All have been declared further down. -->

Note S7: Diagram of all state variables and their dependencies in an exemplar bean simulation
We drew a graph which contains the various state variables in an example simulation and the dependencies among them. Each state variable is simulated by a SimulaObject, here we depicted SimulaConstants, SimulaTables and SimulaStochastic as wedges, whereas all others are depicted as a rounded boxes. The arrows indicate information flow, that is the result of one minimodel goes into the computation of another. The network is strongly dependent on the input file, and somewhat dependent on time, given that computations might switch on given conditions and should thereby be regarded as exemplar. To properly view the graph, enlarge the pdf strongly.