Modelling fungal competition for space:Towards prediction of community dynamics

Filamentous fungi contribute to ecosystem and human-induced processes such as primary production, bioremediation, biogeochemical cycling and biocontrol. Predicting the dynamics of fungal communities can hence improve our forecasts of ecological processes which depend on fungal community structure. In this work, we aimed to develop simple theoretical models of fungal interactions with ordinary and partial differential equations, and to validate model predictions against community dynamics of a three species empirical system. We found that space is an important factor for the prediction of community dynamics, since the performance was poor for models of ordinary differential equations assuming well-mixed nutrient substrate. The models of partial differential equations could satisfactorily predict the dynamics of a single species, but exhibited limitations which prevented the prediction of empirical community dynamics. One such limitation is the arbitrary choice of a threshold local density above which a fungal mycelium is considered present in the model. In conclusion, spatially explicit simulation models, able to incorporate different factors influencing interaction outcomes and hence dynamics, appear as a more promising direction towards prediction of fungal community dynamics.


1.
Introduction. Mycelium, the indeterminate growth form of a filamentous fungus, is comprised of hyphae, i.e. interconnected compartments which create an extending network via their elongation, branching and fusion [21,23,11]. Filamentous fungi interact mainly via competition for space, since nutrients can be acquired from the substratum [1,28]. The main types of fungal interaction outcomes between two heterospecific or somatically incompatible mycelia are the replacement of one mycelium by another (via physical contact, or via chemicals from a distance), and their deadlock, i.e. no mycelium surrenders occupied space to the other [1]. These outcomes of spatial competition are influenced by different biotic and abiotic factors, such as species identity, the presence of grazing inverterbrates, and environmental temperature [1,16,13].
With their interaction outcomes and the consequent community dynamics, filamentous fungi contribute to ecosystem and human-induced processes such as primary production, bioremediation, biogeochemical cycling and biocontrol [14,19,5,20,25]. Incorporating accurate predictions of fungal community dynamics into ecosystem models can hence improve forecasts of such processes, e.g. in decomposition models [18,28]. Despite the literature on the factors influencing fungal community dynamics, and on the importance of these dynamics on ecological processes, only a single study with a simulation model has tested the predictability of empirical community dynamics [17]. Since such simulation models require multiple runs for inference about, e.g., the average behaviour of the modelled system, another attempt could aim for models of ordinary differential equations (ODEs) or partial differential equations (PDEs), which can model the average behaviour of a system deterministically. Moreover, ODEs can model an empirical system more simplistically, for employing the wide range of established ODE techniques. On the contrary, PDE models can provide a more precise description of systems, possibly leading to more accurate prediction of system dynamics. For example, the spatial dependence which spatial PDEs can incorporate could enhance the predictability of dynamics in the complex spatial configurations of mycelia commonly found in empirical communities [8,6].
To our knowledge, only seven theoretical models of interspecific interactions between filamentous fungi have appeared in the literature, a surprisingly small number for a whole kingdom of life [12,9,4,10,3,7,17]. Three of these models require running simulations [12,4,17]. The other four are relatively detailed PDE models, and no prediction of empirical dynamics was attempted by the authors [9,10,3,7]. The first of these PDE models is a system of two PDEs, with one PDE for the concentration of an activator (in the fungal mycelia), and the other for the concentration of a nutrient substrate [9]. The activator diffuses, increases in concentration locally by utilising substrate, but also decays; the substrate diffuses, is captured by the activator, and is also replenished. Without explicitly implementing replacement, this system of two PDEs can not model replacement outcomes, but only 'collisions' of mycelia, i.e. collisions of travelling wave solutions. The second model implements physiological mechanisms, where each mycelium is modelled by five PDEs, with 12 physiological parameters per modelled mycelium [10]. This PDE model is more detailed than the first one of [9], but it can reproduce replacement of fungal mycelia. The third PDE model focuses on only two fungal mycelia interacting in only one spatial dimension [3]. The model is a system of six PDEs, three equations for each mycelium, accounting for: the density of hyphal tips, the mycelium biomass density, and the concentration of an internal growth-mediating substrate. Each mycelium needs eight parameters concerning, for example, hyphal tip speed, anastomosis rate and active translocation cost. Thus, the model is relatively detailed, and limited to only two mycelia on one spatial dimension. Based on the third model [3], the fourth PDE model is also limited to two heterospecific mycelia competing on a single spatial dimension [7].
In the present work, we developed ODE and PDE models which are simpler than the previous models in the literature, and with the aim of predicting fungal community dynamics, validated against dynamics from an empirical system of three species. The basic processes of extension and replacement were incorporated in the derivation of master equations. We took mean-field approximations for the master equations' first moment, to arrive to ODEs for the mean relative abundance of species in well-mixed culture of dispersed mycelia. The non-spatial nature of the ODEs prevented the prediction of even one-species empirical dynamics. The ODE models were then incorporated as reaction terms in reaction-diffusion PDEs. These spatial, PDE models were able to predict one-species empirical dynamics, but suffer from limitations in modelling interspecific fungal interactions, as we discuss in the last Section 4. 2. Empirical system. Three wood decay fungal species were used from the Cardiff University culture collection: Trametes versicolor (abbreviated 'Tv' hereafter; strain TvCCJH1), Vuileminia comedens ('Vc'; strain VcWVJH1) and Hypholoma fasciculare ('Hf'; strain HfDD3). In general terms of competitive ability among these species, Vc is replaced by both Hf and Tv, whereas Hf and Tv cannot replace each other, i.e. they deadlock [14].
Mycelia were inoculated and cultured in 22.4 × 22.4 cm dishes, at 15 • C in the dark, in approximately 0.5 cm thick substratum of 2% (w/v) malt agar. The cylindrical inocula of 0.8 cm diameter were cut from newly developed cultures, and were set with the upper part of the mycelium on the new substratum. The inocula were removed after 2 days (d), to prevent any effects on the dynamics. Reisolation of different mycelial regions showed that visual inspection was adequate for determining mycelial boundaries with unoccupied space and with heterospecifics. Conspecific mycelia were assumed to fuse upon contact of their boundaries [26].
For testing the ability of the theoretical models to predict the dynamics of relative occupancy, the following empirical systems were set. For one-species dynamics, each species was inoculated alone at the centre of a 22.4 × 22.4 cm dish (Fig. 1a), and was allowed to extend unconstrained before reaching the dish walls, or until the whole dish was covered. For three-species dynamics, species were assigned randomly at 49 inoculation sites regularly distributed on the dish (Fig. 1b).
The bottom view of mycelial boundaries was drawn at regular time intervals on transparencies (e.g. Fig. 1b). For measuring species occupancy in time, the transparencies were photographed and processed with ImageJ [24], to obtain the relative occupancy of each species at the time points of measurement (Fig. 1c).
3. Theoretical models. In this Section, simple ODE and PDE models are developed, to predict empirical community dynamics (validated against the empirical dynamics described in Section 2). The models are based on the basic processes of extension and replacement, according to empirical findings on filamentous fungi (Section 1). 3.1. One-species master equation. The empirical system of fungal mycelia on the artificial substrate can be represented as substrate sites occupied or unoccupied by fungus. Thus, the ODE model will describe number or relative abundance of resource sites occupied by the different species. In contrast to the empirical system, this simple ODE model will be non-spatial, assuming each site is neighbour to any other site, such as in a flask with well-mixed, liquid culture of dispersed mycelia [22]. We will assume there is a finite number of available sites, n, for occupancy in the agitated flask. The number of resource sites occupied by species X is symbolised with n X , with unoccupied sites considered an extra species O. For a single species A, the total number n of resource sites in the finite domain of a flask is the sum of the number of sites occupied by species A, n A , plus the number of unoccupied sites, n O : The only basic process at play will be extension of a fungus from occupied to unoccupied resource sites in the well-mixed culture. We will exclude some kind of death or degeneration-recycling processes because we did not observe any in our empirical setting of rich and homogeneously distributed substrate. Species A occupying a site s A can extend to a site s O occupied by (non)species O, resulting in two sites occupied by A: The rate of extension will increase with the relative abundance of unoccupied sites, n O /n, which can be alternatively expressed with Equation (1) as a decrease in the extension rate with the relative abundance of A, n O /n = (n − n A )/n = 1 − n A /n: Hence, the mycelial overall extension rate g A is a function of n A , and e A is the maximum extension rate, when n A = 0. These two extension rates are probabilities of having an extension event per unit of time, i.e. per day in our empirical and modelling settings (unit of measurement: d −1 ).
We will follow an 'occupation numbers point of view' as in [27] for deriving a master equation which is a system of ODEs. Each differential equation concerns the rate of change of the probability Pr(n A ; t) of having n A sites occupied by species A at time t. The master equation will determine the dynamics of the discrete probability distribution for the number of sites occupied by species A. In this onespecies setting, we will have n differential equations, since n A can take the integer values n A = {1, 2, . . . , n}, and we hence have {Pr(1; t), Pr(2; t), . . . , Pr(n; t)}. n A = 0 is omitted because there can be no dynamics without any fungus present initially.
To have n A sites occupied by A at time t + dt, either we had n A − 1 at time t and minimally one of them extended to an unoccupied site during the time interval dt, or we had n A sites at time t and none extended to an unoccupied site. Thus, the probability of having n A at time t + dt is From Equation (2), the probability that one of the n A − 1 sites extends in the time interval dt is g A (n A − 1)dt. Hence, the probability in Equation (3) for all n A − 1 sites is Similarly, the probability that a site occupied by A extends to an unoccupied site in the time interval dt is g A (n A )dt. Hence, the probability of not extending is 1 − g A (n A )dt, and the respective term in Equation (3) for all n A sites is by Maclaurin series expansion. Inserting Equations (4 and 5) into Equation (3), gives Taking Pr(n A ; t) to the left-hand side of Equation (6), dividing by dt, omitting the higher order terms O(dt 2 ), and taking the limit as dt → 0, gives the master equation which can be rewritten, after inserting Equation (2), as

3.2.
One-species ODE model. Since the master equation dictates the dynamics of the discrete probability distribution for n A , we will derive the differential equation for the first moment, i.e. the mean of this distribution. The mean n A (t) of n A (t) in the discrete probability distribution of the master equation is Note that because Equation (9) is a finite sum, Thus, multiplying the master Equation (8) by n A , and taking the sum for all possible n A , leads to a differential equation for its mean: The Pr(n A − 1; t) in the first sum of the right-hand side of Equation (11) can be rewritten as Pr(n A ; t) by changing the limits of the sum from n A = {1, . . . , n} to n A = {0, . . . , n − 1}, and by omitting the n A = 0 term since it is zero. Additionally, the limits of the sum can become n A = {1, . . . , n}, if we add and subtract a term for n A = n. Note that this term is zero because of the n/n division: Doing some calculations in the final sum of Equation (12), gives: The final two sums in Equation (13) can replace the first sum of the right-hand side of Equation (11): Two of the sums in Equation (14) cancel out, leading to Thus, the dynamics of the first moment depend on the second. In long enough time, a fungus will occupy all unoccupied sites in any stochastic realisation, and hence the variance of the stochastic trajectories will tend to zero, Var[n A (t)] = 0.
Since the second moment is n 2 2 . This leads to the final equation for the first moment dynamics: Equation (16) is in terms of mean number of sites in time, n A (t) , and we can cast it in terms of mean relative abundance of species A in time, A(t) = n A (t) /n, leading to the deterministic ODE model of logistic increase in the species mean occupancy (see Fig. 2 for the resulting well-mixed dynamics compared to the spatial dynamics of one-species mycelia extending on empirical dishes): 3.3. Two-species master equation. Although the initially exponential growth of the logistic ODE model does not agree with the quadratic-like growth in the spatial dynamics of one-species mycelia extending on empirical dishes (Fig. 2), we continue towards the development of two-and three-species ODE models, starting from the master equation for two-species. We assume the following basic processes for the two-species system: extension of species A, extension of species B, and replacement of species B by species A (such as Hf or Tv replacing Vc in our empirical system).
As in the one-species case of Equation (2), the extension rates of the two species are functions of the number of sites occupied by A and B, n A and n B : Time ( The rate by which one site occupied by species A replaces a site occupied by B will be a linear function of the relative abundance of B: Again, the extension rates g A , g B , e A and e B , and the replacement rates k A and r A , are probabilities of having an event per unit of time, i.e. per day in our empirical and modelling settings (d −1 ). For the master equation, the probability of having n A and n B at time t + dt will be Pr(n A , n B ; t + dt) · Pr(none of n A replaced one of n B ).
By treating Equations (18)(19)(20)(21) as in the steps with Equations (4-7) for the onespecies master equation, we can derive the master equation for the two species: 3.4. Two-and three-species ODE models. The first moments for the occupancy of the two species are: Following the same approach as in Equations (10-15) for the quantities in Equations (23 and 24), we can derive the differential equations for the first moment dynamics from the master Equation (22): Since species A will eventually replace B in long enough time of any stochastic realisation, the variance of each species will tend to zero, and the covariance of A-B can be neglected. By adopting the mean-field approximations n 2 A (t) = n A (t) 2 , n 2 B (t) = n B (t) 2 , and n A (t)n B (t) = n A (t) n B (t) , Equations (25 and 26) can be written as Equations (27 and 28) are in terms of mean number of sites in time, n A (t) and n B (t) , and we can cast them in terms of mean relative abundance in time, A(t) = n A (t) /n and B(t) = n B (t) /n: For species A, per capita extension is reduced with increasing relative abundance of A or B, and per capita replacement is increased with increased relative abundance of B. For species B, per capita extension is reduced with increasing relative abundance of A or B as well, and its per capita replacement is increased with increased relative abundance of A. From the two-species ODE model of Equations (29 and 30), we can write a threespecies ODE model in which species A and B can replace species C, and species A and B cannot replace each other, as in the empirical system (Section 2): 3.5. One-species PDE model. The one-species ODE model assuming well-mixed conditions could not capture the quadratic-like growth attributed to the effect of space on the dynamics of relative cover or abundance (Fig. 2). Therefore, we can use the one-species ODE model of Equation (17) for a spatial, one-species PDE model in two spatial dimensions, x and y. We cast the model in 2-D since data from the empirical system were of 2-D mycelial boundaries (e.g. Fig. 1b,c). The ODE model of Equation (17) determines the logistic increase in the occupancy of resource sites of a species A. Equivalently in the PDE model, we will assume that the density A(x, y, t) of a mycelium at the x, y point in 2-D space grows logistically in time t (density of species A will be shortly symbolised as A hereafter). A mycelium attains its maximum density, A = 1, with growth rate A (d −1 ). Additionally, the mycelium diffuses with diffusion coefficient δ A (cm 2 d −1 ). This diffusion can be interpreted as the random walk which hyphae follow during extension at the mycelial periphery [2]. The resulting PDE model is the Kolmogorov, Petrovskii, Piskunov and Fisher equation [29]: Since the empirical dish has closed boundaries, we assumed that ∂A/∂x = 0 and ∂A/∂y = 0 orthogonally to the boundaries (Neumann boundary conditions). The modelled dish was a square with sides of 22.4 cm, similarly to the empirical dish dimensions (Fig. 1a).
Starting from an initial density of species A decaying exponentially in the plane, the extending mycelium's boundary attains a constant speed of front c A = 2 √ A δ A [29]. We will throughout fix the ratio A /δ A = 125 cm −2 , to have sufficiently steep initial densities at the inocula and at the propagating fronts [29]. Since the boundary extension rates of species can be known from the empirical system, i.e. c A is known, and assuming A /δ A = 125 cm −2 , we can calculate the parameters A and δ A for a species A. Such parameter estimation appeared satisfactory for one-species dynamics of relative cover with the PDE model, both for the single Tv mycelium (Fig. 3a), and for the three mycelia extending to cover the whole dish (Fig. 3b). Note that to estimate the relative mycelial cover we had to arbitrarily set a density threshold above which a mycelium was considered present locally. A mycelium was assumed present when its density A > 0.01. The relative cover in the PDE solution was estimated by Monte Carlo integration of mycelial presence in space.

3.6.
Three-species PDE model. The PDE model for three species will follow from the three-species ODEs (31-33). Species A and B can replace C with local   Fig. 2a). (b) Three Tv mycelia closely inoculated at the dish centre (inset is a numerical solution of the PDE model at t = 12 d), fusing to form one mycelium which extended and covered the dish (as in Fig. 2b). The PDE model for Tv had initial conditions similar to the empirical setting, with growth rate A = 2.16 d −1 , and diffusion coefficient δ A = 0.017 cm 2 d −1 . It was assumed that a mycelium is present when its density A > 0.01. The relative cover in the PDE solution was estimated by Monte Carlo integration of the mycelial presence (the curve in each panel is 95% confidence region of the mean relative cover in the PDE solution from 100 Monte Carlo integrations).
replacement rates ρ A and ρ B (d −1 ). The rest of the parameters for a species X have been introduced in the one-species PDE model: growth rate of X, X (d −1 ), and diffusion coefficient of X, δ X (cm 2 d −1 ). The resulting model equations are: Again, since the empirical 22.4 × 22.4 cm dish has closed boundaries, we assumed for any species X that ∂X/∂x = 0 and ∂X/∂y = 0 orthogonally to the boundaries (Neumann boundary conditions). The PDE model's dynamics could not agree with the empirical community dynamics (Fig. 4, with Hf represented by species A, Tv by B, and Vc by C). The three-species PDE model's parameter values used to test the predictability of empirical dynamics are given in Fig. 4. The initial conditions for Fig. 4 were exponentially decaying densities of the species at the same locations  Local extension and replacement rates were taken from the mean boundary extension and replacement rates in the empirical dishes, as for the extension in the one-species PDE. It was assumed that a mycelium is locally present when its density is greater than 0.5. The relative cover in the PDE solution was estimated by Monte Carlo integration of the mycelial presence (the curves are 95% confidence regions of the mean relative cover in the PDE solution from 100 Monte Carlo integrations).
as the centers of the inoculated mycelia at the initial state of the empirical dish at measurement time t = 0 d (Fig. 1c).
4. Discussion. The aim of this work was the development of simpler ODE and PDE models, in an attempt to theoretically predict fungal community dynamics under a simple laboratory setting. The ODE model could not predict even one-species dynamics. The one-species PDE model was in good agreement with the empirical dynamics, but its predictive ability for the three-species community dynamics appeared poor. We hereafter discuss this overall poor performance attributed to the underlying assumptions and nature of the models.
We found that space is an important factor in predicting fungal dynamics with our simple models. It is apparent from the logistic form of the ODE (17) that the resulting well-mixed dynamics are different from the spatial dynamics even of a single mycelium extending on an empirical dish (Fig. 2a). The spatial dynamics of empirical relative cover was of quadratic form in the setting of a single mycelium, since we essentially had a disk-like mycelium extending its radius linearly in time. On the contrary, the dynamics in the well-mixed culture of the ODE model tend to be exponential initially, as long as the relative abundance A(t) is relatively small. The same difference in dynamics is apparent and in another example of three Tv mycelia which fuse and cover the whole empirical dish (Fig. 2b). Nevertheless, we went on to develop a two-and three-species ODE model, which helped in the development of the corresponding PDE models implementing spatial processes. Compared to the non-spatial ODE model for one species, the one-species PDE model exhibited good predictive ability against both empirical settings (Fig. 3a,b).
We additionally found that a PDE model of fungal growth and interactions unavoidably suffers from the arbitrary choice of a threshold density which determines mycelial presence. Regarding the initial conditions for the one-species PDE model, all simulated inocula had exponentially decaying density, but we faced an issue when dealing with more than one species. Since in the three-species PDE model we additionally had the replacement rate parameters, we first estimated the extension parameters as in the one-species PDE model. We then assumed that the propagating front of a species X replacing another attains a constant speed c = 2 √ ρ X δ X [29]. The diffusion coefficient δ X was estimated with the extension rates, and we could estimate the replacement rate because the rate of boundary replacement was known from the empirical system. In this way though, the ratio ρ X /δ X < 125 cm −2 , which leads to less steep fronts during replacement. The less steep replacement fronts result in difficulties estimating the relative cover of species, since the arbitrary threshold for mycelial presence must be increased. Thus, the arbitrary selection of the threshold above which a model mycelium is considered locally present, for estimating the model relative cover as in the empirical dish, is an unsatisfactory and likely unavoidable feature of the PDE models. This is particularly significant when the model is used to predict empirical community dynamics which have been recorded in a binary form, i.e. presence or absence of mycelium in each spatial location [17]. If the model mycelia have steep enough boundaries, the prediction of relative cover might be less sensitive to the value of this threshold, although the threshold needed considerable adjustment to attain satisfactory predictions even in the one-species settings.
The PDE approach revealed some additional issues. Setting exactly the same initial conditions as in the empirical setting was difficult because the exponentially decaying densities of the model mycelia could not be initialised exactly, to match the empirical relative cover of the species in each inoculum at t = 0. Moreover, due to the reaction-diffusion nature of the PDEs, it is challenging to measure and parameterise the PDE processes of reaction (mycelial local growth) and diffusion (extension and replacement). Last but not least, there is no way to identify distinct model mycelia of a species with one PDE, and hence it is difficult to extend the present PDE models by relating the extension and replacement rates with, e.g., the individual mycelium cover which has been shown to influence the interaction outcomes and the theoretical predictions of empirical community dynamics [15,17]. An alternative approach would be to model each mycelium with a separate PDE, but this would lead to more involved models, such as previous PDE models on interspecific fungal interactions [9,10,3,7].
In [10], a detailed model with physiology-related parameters has been developed, in which each mycelium is described by a PDE system. It is apparent that this model would be more challenging to parameterise and solve in order to predict empirical community dynamics of multiple mycelia of different species, in comparison to the present PDE model. One way to circumvent these parameterisation and computation challenges was followed in other PDE models [3,7], in which space is reduced to one dimension (1-D), so that mycelia interact on a line. To predict empirical community dynamics, this model must be generalised to 2-D and 3-D, likely making it computationally expensive as well. Lastly, [9] modelled interacting mycelia as local concentrations of activator responsible for converting local concentrations of substrate to biomass with a reaction-diffusion system of PDEs. Due to the nature of this activator-substrate model, there is no way for one mycelium to replace another. Adding appropriate processes for reproducing replacement would lead to a more involved model, limiting the simulation of large enough spatiotemporal scales. Given these limitations, and the issues identified with our simple PDE models, PDE modelling and prediction of empirical community dynamics can be expected to be challenging for realistic systems of multiple mycelia from different species.
In conclusion, ODE models assuming well-mixed resources, and PDE models with their specific features, appeared unable even under our simple empirical setting to predict empirical community dynamics (for a summary of features and performance, see Table 1). Despite their lack of stochasticity and their liability to mathematical analysis, our differential equation models display inherent properties which would not allow them to model and parameterise fungal interactions in space in a plausible way. On the contrary, we believe that spatial simulation models, such as cellular automata or individual-based models [12,4,17], able to incorporate different factors influencing interaction outcomes and hence dynamics, appear as a promising direction for improved prediction of fungal community dynamics [17]. A workaround for differential equation models would be to add more details and levels in an ODE or PDE model, and perhaps this is the reason that previous differential equation models for fungi are more detailed [9,10], or with reduced number of spatial dimensions Table 1. Summary of the features and performance of the models considered in the present work.

Model
Advantages Disadvantages Predictability [3,7]. As a consequence of the details implemented though, such models would become computationally heavier to solve for larger spatial and temporal scales which are characteristic at the level of empirical fungal communities.