Real communities of virtual plants explain biodiversity on just three assumptions

To illuminate mechanisms supporting diversity in plant communities, we construct 2D cellular automata and ‘grow’ virtual plants in real experiments. The plants are 19 different, fully validated functional types drawn from universal adaptive strategy theory. The scale of approach is far beyond that of even the most ambitious investigations in the physical world. By simulating 496 billion plant–environment interactions, we succeed in creating conditions that sustain high diversity realistically and indefinitely. Our simulations manipulate the levels of, and degree of heterogeneity in the supply of, resources, external disturbances and invading propagules. We fail to reproduce this outcome when we adopt the assumptions of unified neutral theory. The 19 functional types in our experiments respond in complete accordance with universal adaptive strategy theory. We find that spatial heterogeneity is a strong contribu-tor to long-term diversity, but temporal heterogeneity is less so. The strongest support of all comes when an incursion of propagules is simulated. We enter caveats and suggest further directions for working with cellular automata in plant science. We conclude that although (i) the differentiation of plant life into distinct functional types, (ii) the presence of environmental heterogeneity and (iii) the opportunity for invasion by propagules can all individually promote plant biodiversity, all three appear to be necessary simultaneously for its long-term maintenance. Though further, and possibly more complex, sets of processes could additionally be involved, we consider it unlikely that any set of conditions more minimal than those described here would be sufficient to deliver the same outcome.


IN TROD UCTION
The continuing level of interest in biodiversity is outstanding (Hooper and Vitousek 1997;Gaston 2000;Tilman 2000;Tilman et al. 2014). References to the term now occur in hundreds of thousands of primary and review articles and tens of millions of web pages.
Because of the direct, indirect and abstract value of plant biodiversity to the human population (Gaston and Spicer 2004;Xu et al. 2020), the search for unifying explanations of it that lie 'beyond description' (Harper 1982) has been intense (Tilman 2000;Silvertown 2004;Gaujour et al. 2012). Some of this theorizing addresses only components of a phenomenon that is undoubtedly a formally complex one (Loreau et al. 2001;Vellend et al. 2013;Wardle 2016). Other approaches, such as unified neutral theory (Hubbell 2001;McGill 2003;Fuentes 2004;Scheffer et al. 2018), are startlingly general and parsimonious.
A difficulty experienced throughout this field is that scale of theory usually outweighs practicality of experimentation, even for very wellresourced research groups. We try here to overcome this difficulty by using well-tried cellular automata (CA) methods (Colasanti et al. 2001;Colasanti et al. 2007;Hunt and Colasanti 2007) in an attempt to find a minimal set of conditions and processes that can support long-term, high biodiversity. We address plant communities alone but, because the methodology is principle-rich, the approach is readily extensible to other areas of biology.
Our approach involves a 'bottom-up' realization (Peck 2004) of different patterns of plant biodiversity in real experimental communities of simulated or 'virtual' plants. We utilize known between-taxon differences in vegetative, physiological and reproductive behaviour as revealed by extensive measurements and observations of plants from a wide range of world habitats. (Grime et al. 1997;Díaz et al. 2004;Pierce et al. 2017). Tests of unified neutral theory are also included.

The hump-backed model
It seems clear that high biodiversity in plant communities involves the mediation, in one form or another, of the so-called hump-backed model (Grime 1973a;Zobel and Partel 2008;Assaf et al. 2011). This model predicts alpha-(local community level) diversity in plant (and animal) communities. It is supported by a wealth of quantitative evidence from field studies involving wide ranges of community productivity within the UK (Grime 1973b(Grime , 1979(Grime , 2001, continental Europe (Kelemen et al. 2013;Cerabolini et al. 2016;Brun et al. 2019), central The postulates of the hump-backed model, as followed here, are due to Grime (1979). Two wide groups of environmental factors are held responsible for its emergence. The first is stress which, in this context, consists of environmental factors that place prior restrictions on plant production, examples being shortages of light, water or mineral nutrients, or non-optimal temperature regimes. The second is disturbance (often a management factor) which causes the partial or total destruction of plant biomass after it has been formed, examples being grazing, trampling, mowing and ploughing, and also extreme climatic events such as wind-damage, frosting, droughting, soil erosion and fire.
The hump-backed model predicts that plant diversity is likely to be low if resources are scarce, if external (destructive) disturbance is high, or if resources are high in the presence of little disturbance. In contrast, diversity is greatest at intermediate combinations of resource and disturbance. These dynamics lead to the model's defining feature: a hump-backed shape when diversity is plotted against either resource or disturbance. In the latter case, the phenomenon reproduces the 'intermediate disturbance hypothesis ' of Connell (1978), (see also Wilkinson 1999;Catford et al. 2012).
In the most general case, a bivariate hump-backed shape emerges when diversity is plotted as a function of total community biomass (Grime 1973a, b;Grime 1979;Rosenzweig and Abramsky 1993;Grime 2001;Pärtel et al. 2007). This relationship goes under several names: diversity may appear as richness, species density or community entropy and biomass as either productivity or production. 'Humpbacked' may also appear as 'humped-back' . We consider hump-backed diversity-biomass to be the most appropriate combination of terms.
This relationship is a very high-level community phenomenon, easily missed (Chase and Leibold 2002;Fukami and Morin 2003;Thompson et al. 2005) by even the largest of field experiments (e.g. Hector et al. 1999). However, it forms only the starting point for the work described here, which comprises very large-scale studies of virtual plants growing in real communities in virtual habitats. By applying combinations of biotic and abiotic environmental conditions to these communities in silico, we test the individual and the combined efficacy of several candidate processes as drivers of the bivariate diversity-biomass relationship.

The C-S-R system of plant functional types
The virtual plants used in our modelling are not representations of any named species or taxonomic groups, they are simplified portrayals known as plant functional types (Shipley 2010;Lavorel et al. 2011;Shipley et al. 2016;Tardieu et al. 2020). Among the many contrasting schemes for the classification of plant functional types, the C-S-R approach taken by Grime (1977) stands out because it delivers considerable powers of prediction from a small number of basic assumptions. These, in turn, lend themselves well to the construction of effective CA rule-bases.
Both the theory and the practice of the C-S-R classification of plant species have been reviewed in detail by Hodgson et al. (1999) and Hunt et al. (2004). Like the hump-backed model itself, the C-S-R system relies upon the combined effects of the two environmental driving factors, stress (principally resource availability) and disturbance. Species are classified into C-S-R types according to their relative success or otherwise when faced with different combinations of these factors. (A detailed slideshow presentation of C-S-R theory and its use within this CA model is given within the Supporting Information.).
Our work makes explicit assumptions which aim to uncover how the hump-backed diversity-biomass shape not only appears to emerge from certain combinations of external driving factors but is mechanistically generated by the interplay of those factors with the functioning of individual plant types. This approach defers to Grace et al. (2014) and Wilson et al. (2019), who called for closer investigation of mechanism in such relationships in preference to a focus on bivariate association.

Our assumptions
2.3.1 First. As a result of past evolutionary trade-offs (Grime et al. 1997;Pachepsky et al. 2001;Díaz et al. 2004;Shipley 2010), plant life has differentiated into distinct functional types (Grime 1979(Grime , 2001Grime and Pierce 2012;Díaz et al. 2016;Garnier et al. 2016;Hodgson et al. 2017). The availability of a pool of different types is assumed to be necessary if plant life is to have the capacity to negotiate the various combinations of resource-stress and environmental disturbance which present themselves within ecological timescales (Chesson and Huntly 1997;Naeem et al. 2016;Körner et al. 2018). 2.3.2 Second. In mixed communities, environmental heterogeneity (Davis et al. 2000;Cannas et al. 2003;Davies et al. 2005;Hodge 2006;Lundholm 2009;Tamme et al. 2010;White et al. 2010;Hillebrand and Kunze 2020) acting at, or immediately above, the scale of the whole plant is an agent that, on both temporal and spatial bases, may retard community processes that would otherwise gain in enthalpy (lose entropy) by moving in the direction of monoculture.

Third. Loss of entropy may also be countered if communities
are open to the incursion of propagules from external sources, such as seed rain or soil seed banks (Thompson 2000;Grime 2001;Gabriel et al. 2005;Aslan et al. 2019).

The general approach
The simulations which we use here ) are experimental systems in their own right (Peck 2004;Hunt and Colasanti 2007;Breckling et al. 2011). Because there is considerable stochasticity in the algorithms, a substantial degree of replication is called for. The simulations are all games played against 'nature' (Lewontin 1961;McNickle and Dybzinski 2013), in the sense that the CA grid simulates externally imposed environmental and management conditions (Ritchie and Olff 1999;Chase and Leibold 2002). The individual 'plants' , ranged across all functional types and present throughout many life cycles, experience a primary filtering by these external conditions and go on to deliver process-based community outcomes (Purvis and Hector 2000;Silvertown 2004;Grime and Pierce 2012). The CA grid also permits plant-mediated environmental factors to operate, which exert a further influence on the emergence of the hump-backed diversity-biomass relationship.
This form of modelling, therefore, does not rely upon abstract associative relations between high-level plant community processes such as diversity, productivity, stability, eutrophication or habitat destruction. It creates bottom-up constructions that are based upon the best available principles that can be assembled from observations of plants at work in the physical world. As the simulations can manipulate many environmental drivers and plant subjects quantitatively and simultaneously, the modelling scenarios created in silico can be colossal in scale.
The CA implementation involves a 2D (ground-plan or chequerboard) model of mixed vegetation, described in detail by Colasanti et al. (2007). A related, more finely scaled, CA model also exists (Colasanti et al. 2001) and this operates mechanistically below the level of the whole plant. It successfully reproduces a wide range of whole-plant-and population-level behaviour . However, the simpler model is the one that better facilitates truly large-scale exploration at the community level. Colasanti et al. (2007) built a model that exploited the properties of the C-S-R system of plant functional types. They drew upon inputs from three sources: Crawley and May (1987), Silvertown et al. (1992) and Colasanti and Grime (1993). Within the C-S-R system, plant attributes ultimately reduce to three broad processes which can, in turn, be used as the functional elements of a CA plant rule-base: (i) the relative ability to harvest resources (vegetative growth), (ii) the relative ability to survive periods of resource deprivation (tissue resilience) and (iii) the relative ability to reproduce by seed (fecundity). These are the operators within the universal adaptive strategy theory of Grime and Pierce (2012), which postulates a fundamental evolutionary tradeoff in the allocation of resources between growth, maintenance and regeneration.

Plant rule-bases
We model 19 functional types here, which together cover the whole range of the C-S-R system in considerable detail (Grime et al. 1988). A probabilistic CA rule-base characterizes the growth, maintenance, and regeneration of each plant type with respect to its sensitivity to variation in resource availability and physical disturbance. All the principles we follow rely upon 'ground truth' supplied by extensive physical observations by Grime et al. (1988Grime et al. ( , 1997, Díaz et al. (2004) and Pierce et al. (2017).
C-S-R space is 3D, but it is uniplanar not orthogonal. The three eponymous types (C, S and R) occupy the apices of this space (see Fig.  1 and Hodgson et al. 1999). In our modelling, the maximal probability of a cell's occupant changing in growth, maintenance or reproduction is P max . For each of those three processes, P max is allocated to the type which is best able to deal with the relevant extreme within resourcedisturbance space. Under conditions of high resource and low disturbance, for example, the competitor (C) is given the P max , for vegetative growth. Similarly, the stress-tolerator (S) is given the P max for survival without access to resource under low disturbance, and the ruderal (R) is given the P max for regeneration from seed under high levels of both resource and disturbance. All the values of P max are numerically the same, whatever the plant process. No plant strategies are viable when low resource availability is combined with high disturbance (see Hodgson et al. 1999 and Supporting Information).
For each of the three allocations of P max , a minimum probability of state change P min is allocated to the type which occupies the position within C-S-R space which is furthest from the type having P max .
P-values for intermediate types are graduated accordingly. It is important that each type's responses sum to unity across all three processes if 'Darwinian demons' are to be avoided and the trade-offs inherent within universal adaptive strategy theory are to be respected. For geometric reasons, this condition can only be met exactly when the P for all types across all processes has a value of 0.3333 recurring, which is the specification of the central type, CSR. However, by adopting a logarithmic gradient between P max and P min it is possible to create the difference P max > 2P min while at the same time relaxing the requirement for summation to unity by less than 3.5 %. Balanced plant rule-bases derived in this way are displayed in Fig. 1. Within each of the three processes, P max takes the value 0.4777 and P min is 0.2325.

Principles of the CA
The modelling uses a 2D (horizontal plan) CA in the form of a 64 × 64 bi-toroidal array of discrete cells (the opposite edges of the square are treated as neighbours in a surface of revolution, thus eliminating any edge effects). Each cell may or may not contain one virtual plant. The whole array is iterated forward in single time steps. The environment that the plant, if present, faces within its own cell at every iteration is described in terms of stress (the probability of the individual being able to access resource during that iteration) and disturbance (the probability of the individual being destroyed during that iteration).
Each cell is a computational object which holds the environmental values which interact with the plant rule-base, and a value that represents its occupancy or otherwise by a one plant type. Values in each cell are updated individually at every iteration and new values are determined from its previous ones and from values within cells in its immediate eight-cell von Neumann neighbourhood.
Within each iteration, every plant is first tested to see whether it maintains its presence within the occupied cell by avoiding a disturbance event. If the plant is not destroyed in that way, then to maintain its occupancy it must either have access to sufficient resources for maintenance or have tissues of sufficient resilience to overcome a temporary 'starvation' .
After all plants have been tested for maintenance, each cell is tested to see if new growth can occur into the cell from plants in any of the eight neighbouring cells. The first part of this test establishes whether there is sufficient resource present for new growth to occur. No growth is possible if this condition is not satisfied. If there is sufficient resource, then growth into the cell may be possible by vegetative spread, local seeding, or incursion from seed rain, depending upon the type of the neighbouring plant. A more generally applied seed rain is also included within treatments described later.

Coding the CA
The cellular automaton is a class-based, object-oriented program implemented in the Python programming language v.3.8.6. The Supporting Information contains a complete and commented listing of the code, including a list of library imports, and the facility for a user to execute specimen runs. The main coding structures which are used to operate the CA are summarized here.
The leading class Cell defines a single cell of the cellular automaton. Each cell possesses attributes in the form of a number of variables: • resource, a float (real) value between 0 and 1.0 that represents the probability of a unit of resource being available; • disturbance, a float value between 0 and 1.0 that represents the probability of the cell undergoing a disturbance event; • xPos, yPos, the x and y coordinates describing the position of the cell within the CA matrix; • occupant, a two-member array which holds an instance of a Two functions operate within each cell at each iteration: • maintain, any plant present is first tested to see if it can survive throughout the current iteration. It risks a disturbance event at a probability determined by the value of the cell's disturbance variable. The higher this value, the higher the probability that the plant will die. If the plant survives the disturbance event, it is then tested to see if there is sufficient resource for it to live on. This is determined by the cell's resource variable, the lower the value the higher the probability that the plant will find it insufficient. However, the plant will only die then if its tissues cannot survive starvation. This is additionally determined by a random function such that the higher the tissue value of the occupant plant, the higher the probability that it will live. The value of tissue is taken from the functional type rule-base ( Fig. 1), which is stored within the program along with the values of all types' grow and seed attributes. • growth, each cell is next tested to see if new growth can occur. The resource variable for the cell, and the grow and seed variables from the plant rule-base, are examined, both for the plant in the subject cell and for each plant in the eight neighbouring cells. If there is insufficient resource, then no growth will occur. The higher the cell's resource value the higher the probability that growth will occur. If growth can occur, each cell is tested to see whether one of three random neighbours can grow into the current cell by vegetative growth. This also involves a random function, such that the higher the neighbouring occupant's grow value, the greater the probability that incursion will occur. The growth from a neighbouring cell will occur even if the cell currently contains another plant or is due to contain one in the next iteration. The selection of three random cells creates an effective balance between vegetative growth and seeding. The number was determined by prior experiment, which found that with fewer than three cells, there was insufficient opportunity for vegetative invasion, and with more than three cells there was too much. If the cell still does not contain a plant in the current or next iteration, the cell is tested to see if one of its eight neighbours can grow into the current cell by seeding. Note that growth by seeding may only occur if the current cell is empty. The higher the neighbouring occupant's seed value, the higher the probability that incursion will occur. The involvement of all eight neighbouring cells in this process was also determined by prior experiment. Finally, if the cell still does not contain a plant in the current or next iteration, there is then the possibility of successful establishment from seed rain. This is also determined by a random function, such that the higher the global value of the variable seedrain, the higher the probability of establishment of a randomly-chosen plant type.
The class Experiment contains the matrix of the cellular automaton. The von Neumann 2D matrix is constructed from pointers emerging from Cell towards its eight Cartesian neighbours. In addition to the global variables resource, disturbance and seedrain, the class involves the following additional variables, which are explained in full detail within the footnotes to Table 1. • timeflexD, the proportion by which disturbance is flexed at each timestep; • timeflexR, the proportion by which resource is flexed at each timestep; • timestep, the interval between instances of temporal flexion; • spaceFlex, the proportional amount that resource and disturbance are flexed across the whole fractal map; • fracd, the fractal dimension (spatial granularity) of the whole matrix.

The base design for all experiments
Nineteen types of virtual plant ( Fig. 1) were established in equal proportions within the 4096-cell array as a spatially random mixture. Some simulations continued for just 50 or 100 iterations (e.g. Figs 2 and 3), but 1000 iterations were more usual. One iteration approximates to the lifespan of the most ephemeral type in the community. Each run began with a new random pairing of resource and disturbance, and this was applied to every cell in the array. One thousand replications were applied to each experiment. Resource and disturbance were chosen within wide probabilistic limits at first, but these were narrowed in some of the later work (see footnotes to Table 1). Even when wide limits were applied, both ranges stopped 20 % short of the lethal extremes of resource = 0 and disturbance = 1. At every iteration, two plant-based outcomes were recorded at the level of the whole array: the Shannon entropy, an index of community diversity (richness and abundance combined, Hill 1973;Spellerberg and Fedor 2003;Jost 2006) and the total 'weighted biomass', calculated as the sum of the number of individuals present within each plant type multiplied by the corresponding 'vegetative' attribute from Fig. 1. When all 19 functional types are equally present, the value of Shannon entropy is 2.94; fewer types and/or unequal presences reduce this value. When one type is present alone, Shannon entropy is zero.
A single experiment with one of these virtual plant communities thus involves the birth, life and death of possibly hundreds of generations of individuals. Each one deals independently with its own resource capture and utilization, and with its response to destructive disturbance, according to its own allocated rule-base.

Additional experimental treatments
Other variables were superimposed selectively onto the base design in order to test the effects of (i) spatial or temporal variability in stress and disturbance, (ii) presence of seed rain, (iii) the assumption of a plant rule-base which reproduces 'neutral theory' . Table 1 gives a schedule of all experiments performed. Footnotes to the table define the way in which each variable was applied and scaled. Executable program code in the Python programming language is available within Supporting Information.

Experiments 1-2: visualization of typical outcomes
Plan views of CA arrays after two specimen runs are shown in Fig.  2. These illustrate a species-rich and a species-poor community. On this occasion, both involve a maximum of just seven functional types. Shannon diversity (entropy) in the two communities is inversely related to total weighted biomass. In Fig. 2A

Experiment 3: the emergence and stability of the diversity-biomass hump
Across a batch of 10 000 runs, each starting with 19 plant types but having a different random combination of global resource and disturbance, the expected hump-backed shape in diversity-biomass emerged after 100 iterations (Fig. 3A). Only 3767 of the runs supported a plant community at this stage. The symbol representing each community is that of the plant type which is most like the net C-S-R composition of the whole (method of Hodgson et al. 1999). Shannon diversity is shown as a function of total weighted biomass in Fig. 3A. Communities resembled type C at the highest total biomasses, where competitive selection operates most strongly, and types S or R at the lowest, where the selection is principally environmental (Wilson et al. 2019). In addition to the hump-backed shape, there were blended transitions between adjacent types across the whole relationship.
The hump-backed shape seen in Fig. 3A is an obvious manifestation of the central limit theorem so, in the absence of more explicit theory, we summarize it using bell-shaped functions such as the Gaussian curve When 'on', all participating plant types followed the non-specialist rule-base of type CSR (see Fig. 1). e In place of spatial uniformity, the array was divided into halves according to a fractally generated pattern. The granularity of the pattern was controlled by fractal dimensions (FD) of −0.6, −0.3, … +0.9 (see Fig. 8). Each replicate run within an experiment began with a new pattern of the same fractal dimension and a new random combination of resource and disturbance. A spatial flex (SF) then moved the values of resource and disturbance randomly and independently upwards or downwards within each half of the pattern by factors of 0.025, 0.05, 0.1 or 0.2. Once set, the FD and SF combination persisted throughout the run. f In place of temporal uniformity, a time step (TS) was introduced by varying resource and disturbance every 2, 4, 8 or 16 iterations. Each run within an experiment began with a different base combination of random resource and disturbance. During a variation, resource and disturbance were flexed (Tf R, Tf D) by adjusting them randomly and independently upwards or downwards from their starting levels by the factors shown. All variations persisted until the occasion of the next variation, which was re-applied with respect to the base combination of resource and disturbance, not the most recent one. g Uniformity of resource and disturbance was applied as in (b), but the probability per iteration that each unoccupied, resourced cell experienced invasion by a single plant propagule from an external source was flexed according to this value. The type of the invading propagule was chosen randomly on each occasion. This curve enjoys the advantage that all three parameters have direct meaning within the relationship being modelled, a being its maximum height, b the location of the maximum within the x-dimension and c a measure of the width of the bell. However, on occasions where the Gaussian curve displays obvious lack of fit, we switch to a four-parameter Rational function: The fourth parameter provides more flexibility to the shape of the bell, if required, although none of the parameter values has any direct meaning within the diversity-biomass relationship. In Fig. 3A, the Gaussian  trendline has a maximum of 1.94 Shannon units at 893 units of community biomass. When the types representing these communities are positioned into the resource-disturbance space from which they emerged (Fig.  3B), they locate themselves in complete accordance with C-S-R theory (see Fig. 1 and the slideshow presentation within Supporting Information). The occupied space is not rigorously triangular, but the horizon between vegetation and non-vegetation is sharp. The plant rule-bases, therefore, capably represent the principles upon which C-S-R theory is founded.
Almost fifty years on, the distribution of types shown in Fig. 3A also exactly reproduces Grime's original formulation of the hump-backed model (Fig. 4, redrawn from Grime 1973a). Diversity was depicted additively then, and types that were later identified separately as S and R were grouped together into a class of species which, while relatively uncompetitive, was nevertheless tolerant of extreme environments or management regimes. Considering these distinctions, the correspondence between Figs 3A and 4 is striking.
Comparing Fig. 3A and B, it is evident that the zone of highest Shannon diversity, i.e. in the very crown of Fig. 3A, corresponds to resource and disturbance combinations that are located just above the centre of the lower margin of Fig. 3B. Within this area if the figure, the net composition of the communities most resembles the types that lie near the centre of the C-S-R classification (Fig. 1). A departure in any direction from this area's combination of resource and disturbance has a markedly deleterious effect on Shannon diversity.
The diversity-biomass hump deflates as the experiment is continued further (Fig. 5A). By 1000 iterations, only 3605 of the original 10 000 communities remain and the hump-backed shape, though still present, has diminished in height. The fitted Shannon maximum is now 0.59 at 885 units of community biomass, so the height of the hump has reduced to 30 % of its former maximum. Most communities have moved towards monoculture (zero Shannon diversity) but, nevertheless, the distribution of representatives within resource-disturbance space remains in accordance with C-S-R theory (Fig. 5B).
As at 100 iterations (Fig. 3B), the same area within resource-disturbance space in Fig. 5B also remains responsible for generating the highest diversity at 1000 iterations. This area is highlighted within Fig.  5B and is enlarged in two stages in Fig. 5C and D.
When environmental conditions are fixed and stable, as in these initial experiments, the combinations which support the highest diversity seem to be very closely defined. In contrast to other regions within resource-disturbance space, which deliver entirely predictable community types (Fig. 5B), the area shown in Fig. 5D is a microcosm supporting almost all shades of community composition. It seems to contain the balancing point of the whole process. In nature, it is perhaps unlikely that a single location could sustain such a precise balance between resource availability and disturbance indefinitely. In the longer-term, therefore, high diversity may involve the intervention of additional factors, such as the heterogeneity and seed rain which are also under investigation here.
The shape of the diversity-biomass relationship remains intact with the passage of time (Fig. 6A) and it also remains centred around 900 units of biomass. Though the earlier scale is quickly lost, a plateau in diversity is approached asymptotically with time and little further decline is evident beyond 600 iterations. These experiments are clearly lengthy enough to approximate to steady-state community outcomes.
The two environmental dimensions, as applied here, run in directions which are numerically opposite in their effect on diversity. As these experiments offer a comprehensive range of environmental combinations, and as many of these combinations are unsuitable to the establishment of vegetation at all (Figs 3B and 5B), it is not surprising to find, from the heatmap shown in Fig. 6B, that at 100 iterations Shannon diversity is maximised along a ridge that runs diagonally across resource-disturbance space. As expected from Figs 3B and 5B, the highest part of this ridge is located at environmental combinations which are low in disturbance and moderate in resource.
The angle which the diagonal ridge makes relative to the resource dimension is steeper than that which it makes relative to disturbance, indicating that (as constructed here) the former dimension is roughly twice as influential in its effect on Shannon diversity. Although diversity diminishes greatly with time, the same route across resource-disturbance space continues to be followed at 1000 iterations ( Fig. 6C). At that stage, the multiple regression which links Shannon diversity (S) to resource (R) and disturbance (D) is (This fit has F = 56.9 and P << 0.00001.) The presence of a negative quadratic term in R indicates that, although the effect of increasing resource is initially positive, diminishing returns eventually result. This is confirmed by the curvature in the low diagonal ridge shown in Fig. 6C. The strands which can be seen emerging from the low peak of diversity in this figure probably arise from performances by individual plant types, which become visible when diversity is very low.

Experiment 4: a test of 'unified neutral theory'
In a repeat of the experiment shown in Fig. 3A (albeit at reduced replication), all plant types were allocated the non-specialist rule-base of type CSR. 'Types' were thus differentiated in name only-a nonfunctional attribute. By 100 iterations, 368 of 1000 original communities had survived but the diversity-biomass shape had become asymptotic, not Gaussian (Fig. 7A). Instances of high total biomass disappeared (cf. Fig. 3A) because the neutral rule-base had deprived the normally most bulky types, in the region of C, SC, CR and their intermediates, of their opportunity to predominate. Across almost all productivities, there was an equal, high diversity among types (names). Only in the most adverse environments (at the left-hand end of the curve) were any types lost during the simulation. Because resource and disturbance were re-randomized in each new run, communities again occupied the entirety of the accessible part of resource-disturbance space. Being functionally indistinguishable, however, the representative 'types' were scattered indiscriminately (Fig. 7B).
Even at 1000 iterations (not shown), no competitive elimination of types occurred at the right-hand side of the diversity-biomass relationship, though losses towards the left-hand side increased and the asymptote lost part of its height. Genuine, between-type functional variations therefore seem to be essential if the diversity-biomass hump is to emerge. Our first main assumption appears to be upheld.

Experiments 5-28: spatial environmental heterogeneity
We superimposed variations in the spatial distribution of resource and disturbance upon the experimental array. Heterogeneity was applied both at different absolute levels and at different levels of fractal granularity (see Fig. 8 and Table 1, footnote e). Throughout, the background levels of resource and disturbance were varied only within their most diagnostic ranges (depicted in Fig. 5C).
Outcomes are expressed in the form of Gaussian curves fitted to diversity-biomass relationships at 1000 iterations (Fig. 9). The object of these experiments was to discover the combinations of variation and granularity that were most effective in raising the height of the Gaussian curve obtained in the absence of spatial heterogeneity (Fig. 5A). Statistical analyses of Gaussian parameter a were used for this determination (Table 2).
It is evident from Fig. 9 that the degrees to which both spatial granularity and resource/disturbance have been flexed clearly succeed in spanning their optimal ranges. Many combinations led to a significant gain in height of the Shannon maximum (P < 0.05), but the two most effective were those having granularity at fractal dimension +0.3 with variation 0.1, and granularity at fractal dimension +0.6 with variation 0.05. These two combinations produced gains of 42 and 55 %, respectively, over the outcome at 1000 iterations in the non-flexed regime, and were taken forward into further experiments.

Experiments 29-64: temporal environmental heterogeneity
To test a second strand of our second main assumption, we introduced dynamic changes into the environmental regime by superimposing variations in the temporal availability of resource and disturbance. Heterogeneity was applied independently to resource and disturbance, both at different absolute levels and at different degrees of temporal granularity (see Table 1, footnote f). Throughout, the background levels of resource and disturbance were again varied only within their most diagnostic ranges (Fig. 5C). Resource was flexed in numerically higher steps than disturbance because of the greater sensitivity of communities to the latter, as already noted. After 1000 iterations and with a time step of eight iterations, the lowest levels of variation in resource and disturbance were the most effective ones at raising the Shannon maximum (Fig. 10). The value 0.70 of achieved by the combination Tf D zero/Tf R 0.0125 represents a recovery of 18 % beyond the position of the untreated community (Fig. 5A). At the shorter time step of four iterations, the most effective variations in resource and disturbance were those of moderate degree. Together with the combination mentioned, Tf D 0.00625/Tf R 0.025 was also taken forward into further experiments.
Taken as a whole, temporal heterogeneity proved to be much less effective in raising Shannon diversity than spatial heterogeneity (Fig.  9). Very few maxima emerged that were significantly above that of the untreated community (Table 3). Time steps of 1 and 32 iterations were also examined but added nothing extra. The prospect of temporal variation alone adding different micro-communities to the array, and hence of raising the diversity of the whole, thus seems to be relatively slight.
Despite versions of both spatial and temporal heterogeneity providing varying degrees of support for diversity at 1000 iterations, it was not possible to approach the height of the full diversity-biomass hump (Fig. 3A) by means of either of these two interventions alone. Heatmap for Shannon diversity within resource-disturbance space at 100 iterations (data from Fig. 3A and B). (C) The same at 1000 iterations (data from Fig. 5A and B).

Experiments 65-73: opportunity for invasion
To test our third assumption (and in addition to the self-seeding possibilities inherent in the base design), we introduced the probability that within a single iteration of the model each unoccupied, resourced cell would experience invasion by a single plant propagule from an external source. The type of each added propagule was determined randomly, since Colasanti et al. (2007) found no appreciable difference in outcome between that form of presentation and one which was weighted according to the fecundity of types already present. Figure 11 shows the effect of an external 'seed rain' (in effect, a composite of both seed rain and soil seed bank) superimposed upon the base experiment at nine different probabilities of invasion. It is inherent in these treatments that high levels of invasion will produce graphical artefacts at low biomass because of the presence of appreciable numbers of recently introduced individuals which never go on to become established members of the community. For this reason, and in the interests of obtaining more informative Gaussian maxima, communities of biomass < 350 units are disregarded whenever a seed rain treatment has been applied.
Despite this, fitted curves at the very highest probabilities of invasion still showed distortion from this cause even with the more flexible Rational function replacing the Gaussian. The spectacular enhancement of diversity that can be seen above level 2E-02 (2 % probability of invasion) should, therefore, be disregarded. The levels 1E-02 and 2E-02 were taken forward into further experiments.
Statistically, all but the very lowest probabilities of invasion produced Gaussian curves which peaked at levels significantly above the collapse depicted in Fig. 5A (see Table 4, P < 0.05). Significant effects began as low as 5E-04, which is one-twentieth of 1 % per iteration, or no better than an even chance of a single propagule of random type being introduced into an unoccupied, resourced cell across all 1000 iterations. This illustrates the high potency of this form of intervention.
Although the height of the full diversity-biomass hump at 100 iterations (Fig. 3A) was approached (Fig. 11, treatment 2.E-02), none of  the levels of application which returned a Gaussian outcome was quite able to equal the former result.

Experiments 74-85: pairwise combinations of treatments
Three groups of pairwise treatments were applied, with two levels (or combinations of level) each in spatial heterogeneity, temporal heterogeneity and seed rain. It is immediately evident from Fig. 12 that the only important interaction between treatments arises when spatial heterogeneity is combined with seed rain. For example, in the case of the leading combination shown in Fig. 12B (the treatment labelled Sp2 + SR2), the Shannon maxima for the two interventions applied singly were 0.92 (Table 2) and 1.55 (Table 4), respectively, but when applied together the maximum was raised to 1.98 (Table 5).
In contrast, no pairwise combination involving a temporal heterogeneity treatment resulted in any multiplicative effect. The outcomes for spatial heterogeneity (Fig. 9) were not changed by adding temporal heterogeneity (Fig. 12A), and neither were the powerful effects of seed rain (Figs 11 and 12C). Davis et al. (2000) postulated a temporal heterogeneity to seed rain interaction, but only the spatial version seems to be effective here.
In Fig. 3A, the Gaussian trendline at 100 iterations has a maximum of 1.94 Shannon units at 893 units of community biomass. In Fig. 12B, the combination Sp2 + SR2, which provides both a spatial Figure 9. Diversity-biomass curves at 1000 iterations under treatments supplying spatial environmental heterogeneity (see Table  1). Different levels of variation (parts A-D) and fractal dimension (FD) are depicted. heterogeneity of 10 % at a fractal granularity of +0.3 and a 2 % probability of incursion by propagules, was able to equal that result by delivering a Shannon maximum of 1.98 (Table 5), at 672 units of community biomass.
Achieving that outcome means that all three of our primary assumptions appear to be upheld (the necessity of functional types, environmental heterogeneity and incursion by propagules). This work has, therefore, achieved its principal aim. Figure 13 provides a full presentation of the result for the combination Sp2 + SR2. Part A shows the outcome at 1000 iterations. The Shannon curve is now almost identical to that of Fig. 3A at 100 iterations, in contrast to the collapsed hump shown at 1000 iterations in Fig. 5A. Admittedly, the left hand-side of Fig. 13A has suffered from the elimination of low-biomass communities, for reasons already described, but the same types of community appear in all areas which are common to Figs 3A and 13A.
With respect to the distribution of communities within resourcedisturbance space, Fig. 13B is comparable to Fig. 5C. Communities of a net composition which resembles type C once again occupy the lower right, R and its close neighbours occupy the upper right, and S and its close neighbours occupy the lower left. This narrow area within resource-disturbance space still appears to contain the pivotal point of the whole C-S-R system, despite the application of the diversityenhancing interventions.

Experiments 86-93: spatial and temporal heterogeneity combined with seed rain
To explore whether further interactions might exist beyond the double combinations, a final set of experiments applied all three forms of intervention simultaneously in combinations of two levels per treatment.
None of the triple combinations added anything to the corresponding double ones because of the weakness of the treatments supplying temporal heterogeneity. Within each of the four sets of pairs of curves shown in Fig. 14, their closeness is due to the two equally ineffectual temporal treatments. The height achieved by each pair of curves (Table  6) is the same as that obtained under their corresponding sets of double combinations ( Fig. 12; Table 5).

Experiments 94-95: two final checks
First, the ultimate stability of the 1000-iteration curve illustrated in Fig. 13A was tested by running 1000 replicates out to 10 000 iterations. The result was virtually identical to that obtained from the 1000-iteration experiment (Fig. 15A) and this was true of all intermediate stages (Fig. 15B). Second, when the 'neutral theory' option (see previously) was superimposed upon the same experiment, the results (not shown) found no interactions between functional type, spatial or temporal heterogeneity, or neutral seed rain.

Caveats
It hardly seems necessary to point out how far removed from physical reality the plants in these experiments are. They are simply 2D 'stickers' which are either present or absent from individual squares on a chequerboard. They contain no explicit representation of physical structure, such as above-or below-ground components, and certainly possess no surrogates for organs, tissues or plant cells. Their morphological, physiological and reproductive properties are represented only with the greatest simplicity, and the different plant types are distinguished only by starkly graduated attribute values.
The environmental conditions with which the virtual plants interact also appear at no lower level than the whole cell of the array and are presented as random locations within just two continua. Between-plant interactions are mediated exclusively by changes in the 3E-02 1.856 n/a n/a 11 4E-02 2.143 n/a n/a 11 5E-02 2.302 n/a n/a Figure 10. Diversity-biomass curves at 1000 iterations under treatments supplying temporal environmental heterogeneity (see Table 1). Different levels of variation (time steps, parts A-D) and temporal frequency are depicted.
environments of neighbouring cells. Community biodiversity is calculated only from a single presence-abundance metric, with community biomass estimated similarly crudely.

Conclusions
But just as these caveats comprise a cost, they also bring benefits in scale. It appears from our modelling that realistic, very long-term patterns of plant biodiversity can indeed be sustained by the presence of distinct functional types and by manipulations in the levels of, and degree of heterogeneity in the supply of, resources, external disturbances and invading propagules. All appear to be jointly necessary and cannot be substituted by the assumptions of unified neutral theory. Other, and possibly more complex, sets of processes might also lead to the long-term maintenance of plant biodiversity, but we consider it unlikely that any combination of conditions more minimal than the set demonstrated here would be sufficient to deliver the same outcome.
Trade-offs in research effort apply just as much as in plant evolution (Harper 1982;Hunt and Doyle 1984). Modelling which aims to reproduce community dynamics over a multiplicity of generations cannot hope also to incorporate all our detailed current understanding of plant morphology, physiology and reproductive biology. To work at scale, 'Simple, high-speed, parsimonious models are required' (Hammer et al. 2019). Plant functional types, for all their approximation and generalization, seem to provide exactly the level of detail that is required to advance towards a mechanistic understanding of previously inaccessible processes at the higher levels of biological organization. Figure 11. Diversity-biomass curves at 1000 iterations under treatments supplying seed rain. The uppermost three curves are Rational (see Table 1). The experiments described here are novel for two reasons. First, they are time courses that simultaneously involve vegetative fecundity, tissue longevity, reproductive fecundity, lifespan, external resource levels, degree of destructive disturbance, spatial and temporal environmental heterogeneity and invasion by propagules. Second, and because of such complexity, the experiments lie utterly beyond the possibility of repetition in the physical world. Typically, our hump-backed curves have arisen from [64 × 64] environmental cells × 1000 environmental combinations × 1000 iterations and thus involve the lives of 4.1 billion individual plants or plant responses. The total of ninety-one 1000-iteration curves specified in Table 1 therefore represents 373 billion individual plant-environment interactions, and the three 10 000-iteration or 10 000-replicate curves add another 123 billion such events. But only at this colossal scale does it become possible to address the entirety of the diversity-biomass relationship.
These experiments illustrate genuine outcomes of the 'ecological theatre and the evolutionary play' (Hutchinson 1965). 'Plant life'  is pitted here against nature. The self-assembling CA communities, though deterministic in detail, are formally complex (Allen and Wyleto 1983) and deliver truly emergent outcomes. As in a game of pinball or roulette, each step towards a final result is entirely lawful, but the result itself is effectively unique and impossible to predict from knowledge available at the outset. Much repetition is necessary to provide clouds of outcomes within which underlying trends can be detected. It is thus not helpful to regard diversity and productivity as processes which are related in their own right (Craven et al. 2020). They are both outcomes of lower-level contests that involve competition and stress/disturbance tolerance. The hump-backed model began life as an associative analysis, then enlarged into a causal, but still fundamentally associative, process (Grace et al. 2014;Gross 2016). A more mechanistic explanation of it is offered here.
C-S-R theory and the hump-backed model are both derived from the same underlying causes, the availability of resource and the prevalence of disturbance. These appear to be the tectonic plates of biodiversity, nothing more than the pre hoc and post hoc conditions within the biotic and abiotic environments which are responsible for selecting successful participants from a pool of available types and determining which, and how many, of these may retain places within the community.
Finally, it is customary in evolutionary thinking (e.g. Denbigh 1989) to regard complexity as a state of high enthalpy (order) and simplicity as a state of high entropy (randomness), possibly on account of the long timescales that are often necessary for the emergence of complexity (e.g. Nerlekar and Veldman 2020). Counter-intuitively, the diversity-biomass relationship does not conform to this pattern (Ray 1994). Low species diversity has the higher enthalpy because it represents an equilibrium within the game against nature in which most plant types have had the opportunity to persist in the community, but in which processes of elimination have worked actively towards a focus on fewer survivors. In contrast, the richer communities are not themselves equilibria displaying high enthalpy: they represent continuing, weak, long-term struggles that never completely resolve. Accordingly, their degree of enthalpy is low.

Future directions in CA
The model is capable of being used predictively in its present form. Because the identities and abundances of plant types within any community may be condensed into a single net location, or 'signature' , within C-S-R space, the distance which separates two communities may be calculated exactly (Hodgson et al. 1999). This could be highly advantageous in investigations (e.g. Van Ruijven and Berendse 2010) which address community resistance and resilience (Nimmo et al.

2015)
. For example, if a stable community were to be exposed to simulated perturbation, the direction and distance through which its signature moves across C-S-R space would be measures of community resistance. Similarly, when perturbation is relaxed the direction and velocity of the signature's movement would be measures of community resilience. Communities inhabiting different parts of C-S-R space could be compared and resistance/resilience mapped across a full range of compositions. 'Drilling down' towards the behaviour of individual types within each zone would then offer testable mechanistic explanations of the two processes. Community succession (Caccianiga et al. 2006;Prach and Pyšek 2009) could also be approached in a similar way.
Extensions of scale and process in plant ecology are also possible by means of CA modelling. CA models at the organ level already exist (e.g. Colasanti et al. 2001;Hunt and Colasanti 2007) and these incorporate far more of plant morphology and physiology than the CA used in the present study. These finer-scale CAs could be extended to accommodate additional trophic levels, including herbivory and predation. Rohde (2005) has examined the utility of CA in ecology more generally.
On the molecular front, fascinating developments of the chequerboard CA could include manipulations of the virtual plants' 'digital genomes', for that is what appears in Fig. 1. Each attribute value there is effectively a 'digital gene' . At present, all genomes are transferred with perfect fidelity into future generations, but this could easily be amended by introducing degrees of imperfection into inheritance, or by effecting transgenesis. Additionally, the regulation of gene expression, currently also perfect, might be flexed probabilistically, or epigenetic effects in the form of acquired variations in gene expression might be introduced.

SUPPORTIN G INFOR M ATION
A PowerPoint explanation of the basic principles of C-S-R plant strategy theory and its implementation in the current model is available in the online version of this article, and an executable Python program showing the development of specimen communities is available at https://colab.research.google.com/ drive/1CKM1rKIsOzdRJsVBaX-GenOr5mQaBWDR.