Since autoregulation of nodulation is essentially a long-distance inter-organ regulatory network, the modelling efforts focus on the organ-scale signalling mechanisms. Functional–structural plant models, which link ‘spatialization of processes in plant functioning and morphogenesis’ (Godin and Sinoquet, 2005), provide an ideal method for this study. With functional–structural modelling, it is possible to simulate the hypothesized signalling mechanisms that are initiated by and affect plant organ development, and then use the regulated plant architecture as a direct reporter to evaluate these hypotheses. At the technical level, there are two major challenges for the application of functional–structural modelling to investigate autoregulation of nodulation: (1) reconstruction of the 3D architecture of legume roots and (2) co-ordination of the multi-rate signalling-development processes.

Unlike the clearly ordered composition and growth of shoots, the root system is difficult to observe and has much more complex patterns. In the past two decades, plant modellers have developed various approaches for collecting root architectural data, handling 3D simulation of root topology and geometry as well as modelling interactions between root growth and environmental factors (Diggle, 1988; Danjon and Reubens, 2008). Most previous efforts focused on woody roots, while legume roots, particularly with nodulation as a distinguishing feature, have not drawn much attention. Since the number of nodules and their distribution are the main phenotypic aspects for studying the underlying regulatory system (Carroll et al., 1985a, b; Delves et al., 1986; Caetano-Anollés and Gresshoff, 1990; Caetano-Anollés and Gresshoff, 1991b; Gresshoff, 2003; Noorden et al., 2006), collecting empirical data on these attributes and reflecting them in the architectural model are crucial to our work. Non-destructive or automatic technologies for root data collection, including those using Ground Penetrating Radar, CT imaging and 3D laser scanning (Danjon and Reubens, 2008), still have limitations in resolution and branch detection or have restrictions due to environmental conditions, which lowers their effectiveness or feasibility in application to this study. The semi-automatic scanning method described by Lira and Smith (2000) could help in counting nodules in high-resolution root images and could possibly be improved to recognize nodules automatically. But the reliance on 2D images taken from one angle means that the 3D lateral branching cannot be well-classified by this technology and some nodules hidden by primary or lateral roots would be ignored. On the other hand, 3D digitizing approaches (Room et al., 1996; Sinoquet and Rivet, 1997) could position plant organs precisely. However, legume roots are usually highly flexible and therefore change their 3D position once removed from the soil, which makes it difficult for digitization to capture their spatial patterns.

Recreating the architecture is not enough; an architectural model is more useful when it helps to reveal internal and environmental factors that influence its development, which requires the linking of plant structure and function. Most previous functional–structural models of root development (Danjon and Reubens, 2008) take environmental factors into account. However, the root architecture is ‘shaped’ not only by environmental factors but also by the influence of endogenous signals such as hormonal stimuli (Aloni et al., 2006). For functional–structural modelling of internal signalling control, the stage has been set by previous work (Janssen and Lindenmayer, 1987; Prusinkiewicz and Lindenmayer, 1990; Buck-Sorlin et al., 2005, 2008; Lucas et al., 2008) using discrete information transfer based on context-sensitive L-systems. This work is built on through the integration of internal signal regulation with root architectural details, the visualization of signal allocation in the root system and, in particular, the synchronization of signalling and developmental processes with various empirical or hypothetical rates.

To characterize lateral roots based on this emission pattern and collect their developmental data, a ‘RULD’ root mapping method (Han et al., 2007) was developed. With the ‘RULD’ method, the first-order laterals are categorized into quadrants (right, R; up, U; left, L; down, D), according to the relative positions of their emission points to the obtuse angle composed by the two cotyledons (Fig. 2). Then the lateral roots are further classified and identified by ‘regions’ (a 50-mm-long section on the primary root), ‘segments’ and ‘sites’ (Fig. 3).

To capture nodule distribution information, a ‘nodulation section’ (Han et al., 2009) is defined for each region of the primary root and for each first-order lateral root (Fig. 4). The nodulation section covers the distance between the first and the last nodules in a particular primary-root region or on a lateral root. The relative locations of the first and the last nodules to the starting point of the primary or the lateral root are measured to position each nodulation section. The length of a nodulation section and its number of nodules are used to calculate the nodule density for this section.

The empirical architectural data collected with the above measurement methods are then used to drive the root simulation. The primary root architecture is extended by the root tip, where potential positions for nodulation or lateral emission are also created. The root elongation has an appropriate rate to match the empirical length value, which defines the corresponding root tip location. There are three major rules for modelling the primary root elongation. The potential sites pre-set by rules (1) and (2) are only spatial markers; whether and when nodules or lateral roots will be formed there depend on other conditions derived from empirical data. The maximum number of nodules formed in each primary root region during each day is used to restrict further nodulation. At a potential nodulation site, a nodule will be formed if, and only if, the relevant maximum number of nodules is not exceeded. Whether a lateral root will be emitted from a potential formation site is determined by branching probabilities, while how soon the lateral will appear and how far it can elongate are controlled by the empirical daily growth rate. The apical zone around the root tip in the architectural model has no nodules attached and remains unbranched, as the maximum nodule number and the lateral growth rate for this area are both 0.

The architectural model was built with an L-system-based language ‘cpfg’ in L-studio (Prusinkiewicz and Lindenmayer, 1990; Hanan, 1997; Prusinkiewicz, 2004). In L-system plant models, a shoot organ is usually represented by a single module that not only contains its developmental information but also simulates its structure graphically. A number of such modules are then assembled into a string to represent the entire shoot architecture composed of different organs. The dynamics of plant development is supported by step-by-step application of a set of rules called ‘productions’. The productions are usually written as where ‘predecessor’ is an original module while ‘successor’ is a module or a string of modules to replace the predecessor as long as the ‘condition’ is matched. At each simulation step, the modules in the current string are matched against the productions and only those with matched predecessor and condition are applied. If multiple productions have the same predecessor, they are checked sequentially and only the first one with matched condition is allowed to produce its successor. After all modules have been processed, the resulting string represents the structure of the plant at the end of the time step.

In this root study, the elongation information (such as overall root length already created) is recorded within a module representing the root tip, while the graphical role is played by a set of ‘sub-modules’ each with the same length (UNIT). The standard sub-modules here are similar to the elementary units used to model oil-palm root system (Jourdan and Rey, 1997). During root elongation, a sub-module is added behind the root tip module at each simulation step, forming a string of sub-modules to make up root structure. For example, if ‘RT’ is used to represent root tip, ‘R’ to represent a sub-module, the initial structural string is ‘RT’ and the production is then the string will become ‘R RT’ after one step and ‘R R RT’ after two steps. The sub-module length – UNIT – is definable by the modellers or users depending on specific requirements in specific cases. In this study, to allow the increment of root length to match conditions for setting potential nodulation or lateral formation sites, the value of UNIT should be smaller than the minimum interval between two successive nodules or laterals. For details of the time scale, see the next section. The algorithm for primary root elongation with sub-modules is illustrated in Fig. 5. Its implementation using the ‘cpfg’ language is given in the text file in Supplementary data (available online). Although all the sub-modules have the same length, their widths vary and increase with radial growth, capturing appropriate diameters at different stages of root development. The lateral root elongation is simulated in a similar way except that the second-order laterals are not taken into account. The simulation of root heading behaviours is supported by methods from the ROOTMAP model (Diggle, 1988) and is representative only. Since the plants used for studying autoregulation of nodulation in this research are at their early developmental stages, the mortality of lateral root tips is not considered. A sample visualization of soybean root architecture with nodulation is given in Fig. 6 and the video in Supplementary data (online).

To evaluate the root reconstruction, growth data of a soybean hypernodulation genotype nts1116 (Carroll et al., 1985a; mutated at V837 of the GmNARK receptor kinase gene) are incorporated into the architectural model. nts1116 was shown to have severely reduced in vitro kinase activity (Miyahara et al., 2008), consistent with a 3x elevated nodule number. The nts1116 plants were cultured in glasshouse conditions over a 16-d period, inoculated on the second day, with five plants sampled for destructive measurements every 2 d starting on the third day. The comparison of real-plant versus virtual-plant root systems, on lateral root branching and nodule distribution (Han et al., 2009) as well as on root length, indicates a good fit between empirical data and simulation results (Fig. 7).

Since there is more than one signal involved in autoregulation of nodulation and the signalling and developmental processes have various rates, a co-ordination mechanism to synchronize these multi-rate processes has been developed. Firstly, each signalling or developmental event has two states: ‘activated’ and ‘stopped’. During a single L-system time step, only one sub-module can be added to the current structure or be passed through for signal transport if the relevant elongation or signalling event is activated. When their states are marked as ‘stopped’, those signalling or developmental events do not occur and will wait to become activated, independently. The difficulty is in switching between these two states for each event so that the signalling-development processes can be dispatched in a synchronized way. To achieve this, a concept of ‘time division’ that divides one day's time into lower-scale time sections (e.g. one day's time can be divided into 24 h – each hour is represented by such a division) can be defined. Each signalling or developmental event e_i is assigned a certain number of time steps, defined as c_i, during each time division. Assuming the rate of e_i is r_i and the number of time divisions per day is DIV, the value of c_i can be calculated using eqn (1): where ceil(x) is a function returning the smallest integer value not less than x. The number of time steps allocated to the fastest signalling or developmental event during a division, equivalent to max{c_i|i = 1, 2, … , n}, is the total number of time steps this division consists of. That is to say, the event with the fastest rate keeps occurring over all time steps during the division, while the events with slower rates are only activated over a smaller number of steps and are ‘stopped’ during the remaining steps of a division. When a new division is initialized, all the stopped events will become activated again. Therefore all processes are synchronized at the beginning of the division. The synchronization process is illustrated with two sample events in Fig. 9. This strategy allows as many signalling-development events to be involved and synchronized as possible.

With the integrated and synchronized signalling and developmental processes, the functional–structural model of autoregulation of nodulation can produce different signal allocation and regulation patterns based on different parameter settings (Fig. 10). For example, by colouring each sub-module according to its SDI concentration value, the allocation of SDI at different parts of the root can be visualized in detail (Fig. 10A–C), while the inhibited nodules (that are not apparent in nature) can either be visualized (Fig. 10D–F) or filtered (Fig. 10G–I) to help in analysing the nodulation pattern. These functionalities could be used as the basis of virtual experiments to investigate unknown or unclear attributes of autoregulation of nodulation.

Yet for legume plants, the root architecture is composed not only of primary and lateral roots, but also of nodules. The nodules play a critical role in supporting legume plants' growth and in providing fixed-nitrogen to the ecosystem, so must not be ignored. In this paper, the current methods to capture nodulation patterns from the plants are presented and recreated in computational models.

The reconstruction techniques developed in the present study have some limitations. The data collection and model evaluation are based on young soybean roots, where the second-order laterals and mortality of root tips are not considered. These aspects of root reconstruction will need to be addressed if the study is expanded to older plants grown in field conditions. At this stage, the modelling focus has been on soybean. For other legume roots, the primary root elongation algorithm is also applicable, but modification of the data collection method may be required depending on whether the first-order laterals there can be characterized into quadrants (or whether the radial angles of lateral emission are around 90°).

The observable plant structure serves to report the unobservable regulation mechanisms. Using functional–structural modelling, long-distance signalling for autoregulation of nodulation has been integrated with architectural development. To effectively co-ordinate the signalling and developmental processes with various rates, a synchronization algorithm based on the use of sub-module structure and context-sensitive L-systems has also been developed.

The functional–structural model, integrating root architectural details with signalling control, enables parameterization of signalling hypotheses and produces regulation patterns within different forms (e.g. signal allocation and nodule distribution). This provides a basis for implementation of virtual experiments to analyse or evaluate mechanisms of autoregulation of nodulation that are unclear.