Team:DTU-Denmark/GrowthModelling
Growth Modeling
Under the lens
Filamentous fungi have a very characteristic growth pattern on solid medium, forming a large network of interwoven filaments (3). If we can accurately predict how a fungus will grow over time, it will allow us to determine the optimal growth conditions to gain the best end product. Therefore, we have investigated how the fungal morphology changes over time at two levels of details.
Animation of fungal hyphal growth. If you are interested in a more detailed explanation, click here. 
Microscopic view
One of our models focuses on the morphology at the hyphal level by simulating the movement of hyphal tips, branching rates, extension rates and the density levels during a growth period in two dimensions. All of the code scripts can be downloaded in this zip file.
The mycelium of a fungus consists of many interwoven hyphae, and the density depends on how many filaments there are in a given location. In Fig 1 below, it is possible to see microscopic pictures of fungal mycelium. Three different levels of zoom illustrate how the network looks, where it can be observed how they interlink and how a fungal filament can branch into more.
Fig. 1: Snapshots of mycelium development of Aspergillus oryzae. These are representative microscopic images of how a network of intertwined hyphal filaments could look in a microscope.
Simulation of the mycelium development
Fungal growth is initiated by $n$ number of spores, and a branch will start to extend from each of the spores added to the space. Following along one of these branches originating from a single spore, the hyphae will grow in a direction $\theta$ with a tip extension $r_{tip, i}$. A branching event, in which a new branch is formed from the first branch, can occur with a probability $q$.
Tip extension rate is calculated by using the equation below, which considers growth kinetics for the fungi and the amount of substrate available. It essentially outputs the accelerated growth dependent on the amount of substrate available, where the accelerated growth equation depends on fungal kinetics and branch lengths (1, 3).
\begin{equation}
r_{tip, i} = \bigg(k_{tip,1} + k_{tip,2}\cdot \frac{l_{br,i}}{l_{br, i} + K_t}\bigg)\cdot\bigg(\frac{S}{S+K_s}\bigg)
\end{equation}
$k_{tip,1}$ is the initial tip extension rate of the branch and $k_{tip,2}$ is the difference between the maximum extension rate and $k_{tip,1}$ (1, 3). $K_t$ symbolizes the time it takes to reach half of the maximum extension rate. The length of branch $i$ described by $l_{br,i}$. $S$ is the substrate concentration and $K_S$ corresponds to the substrate concentration to reach half of the maximum growth level (4).
As the starting coordinates of this simulation are $(x_0, y_0)$ and end coordinates $(x, y)$, the length of branch $i$ can be calculated as the distance between two points:
\begin{equation}
l_{br, i} = \sqrt{(x-x_0)^2 + (y-y_0)^2}
\end{equation}
By dividing the growth area into a grid of $w\cdot h$ areas, it is possible to investigate the uptake of substrate, hyphal movement and the development of biomass through the simulation.
At the start of the simulation, it can be assumed that the initial substrate concentration $S_0$ will be distributed evenly across the grid. When the spores are placed randomly, the branches start to develop by using the substrate available around the hyphal tip. So for each branch, it is checked which tip is in which grid and whether there is any substrate available. That means that one can track the substrate depletion in each grid, and when there is no more substrate in one of the areas of the grid, the hyphal tips located there can no longer grow and can be considered as inactive hyphae. If the total area around the mycelium no longer contains substrate, the growth will simply stop.
In the same way as tracking substrate depletion, biomass production can also be viewed in a grid. But instead of calculating a loss of substrate, each new branch extension is considered a gain in biomass or the density $d$. It should show the same mechanics as the substrate depletion, as these two are directly related to each other by the equation $l_{br,i}$.
As the fungi grow in size, the computational cost also increases. We were limited in running the simulation due to the power of our computers, thus resulting in us introducing two restrictions: the number of hyphae $M$ and the number of steps in the simulation $N$.
Results
The simulation ran with the parameters listed in table 1, which painted a detailed picture of the time history of growth from the spore to the densely branched network of hyphal filaments. The parameters are based on growth of Aspergillus oryzae in Spohr (2).
| Parameter | Value | Source or rationale |
|---|---|---|
| $k_{tip,1}$ | 80 $\mu m\cdot tip^{-1}\cdot h^{-1}$ | Taken from Spohr (2) |
| $k_{tip,2}$ | 75 $\mu m\cdot tip^{-1}\cdot h^{-1}$ | Taken from Spohr (2) |
| $K_t$ | 5 $\mu m$ | Estimated |
| $S_0$ | 50000 mg/L | Estimated |
| $K_S$ | 200 mg | Estimated |
| $M$ | 100000 | Set as a simulation limit |
| $N$ | 5000 | Set as a simulation limit |
| $q$ | 00.05 | Estimated |
Table 1: - Parameters in a typical simulation run.
The end results can be viewed in Fig 2 below, where three animations are shown of the hyphal development: hyphal movement and locations, biomass development and substrate depletion. It is the same mycelium simulation in all three animations, where we can observe that the density increases as the substrate level decreases.
Fig. 2 Left: Hyphal development over time from 10 spores. Middle: - Density development as the hyphae grows. The darker the color the higher the density. Right: - The depletion of substrate over time, as the fungus uses the substrate to extend its hyphal network. The lighter color indicates lower levels of substrate.
Fig 3 shows how the hyphal lengths are increasing in the beginning of the growth simulation, whereas there are initially few hyphal tips. But as time progresses, the number of hyphal tips starts to increase drastically due to there being more and more hyphae and the original branches being longer. The frequencies of hyphal lengths and branching events can be seen in Fig 4, where the majority of the hyphal elements are very short and haven’t branched that much. Those that are longer are also those which have branched more often than the rest. The substrate level decreases as the amount of biomass increases, where this pattern can be observed in Fig 5. All of the substrate is not gone, but that is due to the fact that the mycelium growth was not that efficient in using all the substrate in the given simulation time.
Fig. 3: The number of hyphal tips (red) and total length of all hyphae (black) over time in the simulation.
Fig. 4a: The frequencies of hyphal lengths, where it can be seen that the majority of the lengths are lower than $\mu$m.
Fig. 4b: The frequencies of the number of branching events per hyphae, where most of the branches only have less than 5 events and those that have more are the oldest branches in the simulation.
hi
Fig. 5b: The increase in mycelium biomass over time, where the curve follows an exponential growth mechanism.
Fig. 5c: The substrate level over time decreases as the mycelium develops.
hi
Substrate dependent growth
The fungus can grow and fill out the space since the substrate concentration is sufficient to support this pattern. If there is not enough substrate available, it is impossible for the fungus to fill out the space. This can be illustrated in the simulation if there is not enough substrate available. For example in fig. 3, the substrate concentration is a fifth of the level in fig. 6. As can be observed, most of the hyphae uses up all the substrate available in their grid area, meaning they do not grow into another grid area as they are locked in.
Fig. 6: An example of substrate limiting growth, where the hyphae uses all available substrate in their squares and are therefore unable to escape to areas with more substrate.
Future development
The model shows how the fungus grows in 2D, so it can be interpreted to work as looking at a petri dish in detail. If the model were expanded into 3D, the result could be related to filling out a form and still studying the hyphal interactions at a very detailed level. 3D models of fungal growth do already exist, for instance see Lejeune (3) that simulates the growth of the filamentous fungi Trichoderma reesei.
There are many different factors that influence the growth kinetics and this model only includes a few of them. Parameters such as temperature or oxygen levels could be implemented to get the simulation to work more realistically.
A Simple PDE model of Fungal Growth
We now wish to move from a model describing the individual hyphal growth to a continuous model, which can describe the change of biomass density as a function of time and space. This can be seen as a meso level model, which only needs information about the substrate level, conversion efficiency and internal diffusion rates to model the spread of the fungi through time and space.
First, the model will be introduced, along with a table of relevant parameters. Afterwards some key results will be presented.
The Model
The model presented below is based on Falconer (5), but greatly simplified because of the non-complex morphological behaviour of our fungi. It should be noted that the full Falconer model was implemented, but all relevant effects can be recreated with the simplified model. \begin{equation} \frac{\partial b_i}{\partial t} = \zeta D_b \nabla^2 b_n \end{equation} \begin{equation} \frac{\partial b_n}{\partial t} = (1-\zeta) D_b \nabla^2 b_n \end{equation} \begin{equation} \frac{\partial n}{\partial t} = D_n \nabla^2 n +s(\lambda_1b_n+\lambda_2b_i) \end{equation} \begin{equation} \frac{\partial s}{\partial t} =-(\lambda_1 b_n+\lambda_2 b_i)s \end{equation} The four equations describe how the fungal biomass and the substrate change over time. The fungi is modeled as consisting of mobile and immobile biomass, the immobile biomass can be either non-insular ($b_n$), corresponding to the hyphal tips of the colony that are capable of large scale extraction of external resources. The mobile biomass ($n$) represents the uptake of nutrients from the available substrate, and is used to regulate the diffusion rate $D_n$, which is lowered to $10^{-7}D_b$ if $n>n_0$, ie. the biomass exceeds a critical value: The parameter $\zeta$ is the conversion rate of non-insulated biomass to insulated biomass. The uptake efficiency of the mobile and immobile biomass are denoted respectively $\lambda_1$ and $\lambda_2$.
| Parameter | Value |
|---|---|
| $\lambda_1$ | 0.95 |
| $\lambda_2$ | 0.01 |
| $D_b$ | 10 |
| $D_n$ | 10 |
| $\zeta$ | 0.01 |
| $n_0$ | 1 |
The system was spatially discretized on a 150x150 grid, which created a set of ordinary differential equations for each grid point (for a more thorough introduction see here). These were numerically integrated using a numerical $2^{nd}$ order strong stability preserving Runge Kutta method, using the package DifferentialEquations.jl in Julia (6). No flux boundary conditions were imposed, which is also known as a Neumann boundary condition, these are deemed realistic since we will be growing our fungi in molds. Alternatively periodic boundary conditions could have been interesting if it was decided to make a spatial grid of inocolum. All calculations were vectorized and parallelized, which allowed us to simulate large systems on personal laptops.
Simple Radial Growth
Loading the density obtained from the micro model in a 9x9 grid, and using that as our initial mobile biomass, we can see how the growth of the initial spores evolves in the meso model. The initial mobile biomass was placed centrally in a 150x150 times grid, and an initial homogeneous concentration of substrate was used.
Fig. 7: An example of the most simple fungal growth, with an initial central growth which spreads outwards radially. 
Note how the density maxes out, while the fungus expands radially outwards, akin to what is observed in a lab setting.
Substrate Limited Growth
As described in the micro model, the meso model of course also limits the growth of the fungus based on the available substrate. This naturally occurs since any increase in fungal biomass is caused by a conversion of substrate to biomass. The model predicts a growth similar to real measured growth of fungi (7), although the scaling of time and biomass are system dependent the overall trend is the same. These results are similar to those obtained in the micro model, and shows that the model can also simulate realistic growth curves in the simple case. Note the asymmetry of the final growth, caused by an initial asymmetry of the initial spores grown in the micro model.
Fig. 8: The fungal biomass is shown on the left after 100 timesteps, and on the right the total fungal biomass is shown as a function of the time. Note the lag phase, going into the exponential phase, corresponding to rapid expansion, and lastly settling in and stopping the rapid expansion because of the limited space and nutrients. 
Several Colonies
In the final product, it will be interesting to investigate how several different colonies can fill up a space, and investigate how the distribution of fungal biomass will be in the product. This model allows us to start several different colonies, by setting the initial concentration of mobile biomass as non-zero in several places and letting the system evolve.
Fig. 1: No antagonistic behaviour between colonies has been observed in S. commune
Fig. 1: A system with two colonies starting out next to each other and merging.
The two fungi met after a few timesteps and grow together to almost form one. As the substrate concentration is lower on the outskirts of the fungi, the periphery of the colonies have a lower density than the center.
This model only works for non-antagonistic colonies (8), since otherwise, we would need a growth inhibitor when the two colonies meet. But as the picture of S. commune show, and based on observations of the growth of A. oryzae, we did not observe antagonistic behavior between colonies of the fungi we looked at.
Future Work
Expanding the model to 3-spatial dimensions should be fairly trivial, and would offer a better representation of how the density of mycelium changes as a function of time in the final molds. But the primary challenge is verifying that the density distribution is as predicted by the model, this could perhaps be achieved by measuring absorbance of an agar plate with fungal growth. Furthermore the parameters of the model have not been experimentally measured.
Conclusion
The micro model is already based on parameters from the growth rate of A. oryzae, but could potentially be used for other species of fungi. None of our experiments measures the branching rates and how the hyphae interact, but the model would be more precise if these were included. The model itself allows us to gain a good understanding of how the fungal behaves at very detailed levels.
The meso model could be calibrated based on the measured growth rate of whatever fungi is used, as long as there are no antagonistic behaviors between colonies, and an arbitrary shape could be filled to a certain minimum or average density. Since our statistical model shows that the density of fungi is an important system parameter, this allows us to estimate the fungal density of our boxes and to calculate how fast different shapes are filled with mycelium.
(1) Lejeune R, Nielsen J, Baron GV. 1995. Morphology of Trichoderma reesei QM 9414 in submerged sultures. Biotechnol Bioeng 47:609-615.
(2) Spohr A, Dam-Mikkelsen C, Carlsen M, Nielsen J, Villadsen J. 1998. On-line study of fungal morphology during submerged growth in a small flow-through cell. Biotechnol Bioeng 58:541–553.
(3) Lejeune R, Baron GV. 1996. Simulation of growth of a filamentous fungus in 3 dimensions. Biotechnol bioeng 53:139–50.
(4) Monod J. 1949. The Growth of Bacterial Cultures. Annual Review of Microbiology 3:371–394.
(5): Falconer RE, Bown JL, White NA, Crawford JW. 2005. Biomass recycling and the origin of phenotype in fungal mycelia. Proceedings of the Royal Society of London B: Biological Sciences 272:1727-1734.
(6): Rackauckas C, Nie Q. 2017. DifferentialEquations.jl – A Performant and Feature-Rich Ecosystem for Solving Differential Equations in Julia. Journal of Open Research Software 5:15.
(7): Held P. 2010. Monitoring growth of beer brewing strains of Saccharomyces cerevisiae. Biotek Application Note:1-6.
(8): Falconer RE, Bown JL, White NA, Crawford JW. 2008. Modelling interactions in fungi. Journal of the Royal Society Interface 5:603-615.