Prediction and experimental evidence of the optimisation of the angular branching process in the thallus growth of Podospora anserina

Based upon apical growth and hyphal branching, the two main processes that drive the growth pattern of a fungal network, we propose here a two-dimensions simulation based on a binary-tree modelling allowing us to extract the main characteristics of a generic thallus growth. In particular, we showed that, in a homogeneous environment, the fungal growth can be optimized for exploration and exploitation of its surroundings with a specific angular distribution of apical branching. Two complementary methods of extracting angle values have been used to confront the result of the simulation with experimental data obtained from the thallus growth of the saprophytic filamentous fungus Podospora anserina. Finally, we propose here a validated model that, while being computationally low-cost, is powerful enough to test quickly multiple conditions and constraints. It will allow in future works to deepen the characterization of the growth dynamic of fungal network, in addition to laboratory experiments, that could be sometimes expensive, tedious or of limited scope.

Filamentous fungi are characterized by their ability to form an interconnected hyphal network, the mycelium, based upon some key cellular processes, i.e. hyphal tip (or apex) growth, branching and hyphal fusion (also known as anastomosis), thus conferring a great flexible morphology and a remarkable capacity of adaptation to very diverse ecosystems 1 . In particular, saprophytic fungi, known as short-range foragers, evolve in a highly competitive habitat and thus appear to be highly challenged by resource-limited and patchy environment, combined to a fierce competition with other organisms 2 . Then, these fungi must adapt their growth and have always to find a compromise between the need to occupy the space potentially threatened by other organisms (colonization) and the need to optimally draw resources from where the thallus has developed (densification). Previous studies have reported that hyphal tips at the biomass edge could be associated with exploration while hyphal tips behind the biomass front could be most associated with resource exploitation 3 . In the same way, Lew 4 described the apical cells (or leading hyphae) from the fungal network as the first cells to invade new territory and are generally engaged in nutrient acquisition and sensing of the local environment. However, until now, the way a fungus optimizes its growth through an efficient compromise between the maximization of the surface occupancy and the increasing production of length is still unclear.
For many years-a review can be found in 5 -mathematical modelling and measurement of specific quantitative parameters through the use of image analysis have greatly contributed to a better understanding of fungal network expansion and topology. Usually, these approaches attempt to explore the fungal development in two or three dimensions, in order to define some specific macroscopic and microscopic observable and to describe, explain and predict development of the fungal network in various more or less constraining environments. It was shown in 3 that, in the saprophytic fungus Rhizoctonia solani, translocation (nutrient transport) was mainly diffusive in homogeneous environment, when growth experiments were compared to a detailed mathematical modell; in the latter case, translocation was considered to have both diffusive and metabolically-driven components. Falconer et al. 6 demonstrated that fungal phenotype could be modelled as an emergent phenomenon resulting from the interplay between local processes involved in nutrient uptake and remobilization of internal resources, and macroscopic processes associated with their transport. In order to enlarge the spatio-temporal zation of branching network 13,14 : (1) Random networks for which the connections between the nodes are random and are free from the physical distance that separates them; (2) Totally interconnected networks which allow all the nodes to be linked together; (3) Freescale networks which are characterized by a number of connections between nodes which follows a power law (i.e. many nodes which have few connections but some nodes which have many); (4) Smallworld type networks which are characterized by strongly interconnected subnets, which are themselves connected to each other via a few privileged links. It is also necessary to distinguish between highly centralized networks, which assume privileged nodes in order to transmit information from one node to another (the geodesic in terms of information transmission is then not optimal), and networks without privileged centers that allow an optimal transmission of information. A third important aspect is the existence of networks that do not have spatial dimensions (internet for example) and others whose spatial dimension is critical in the sense that the probability of connection between two spatially distant nodes is zero. Finally, we must distinguish between static networks (in the sense of the number of nodes and not of the activity that takes place there) and dynamic networks which are capable of complicating and increasing their number of nodes and links as a function of time.
The study of P. anserina in 12 allowed to identify the core characteristics the simulation should be based on. P. anserina thallus is a spatialized and dynamic network (Fig. 1). Accordingly, the probability of connecting spatially distant nodes (or vertex) is zero. However, if the thallus of P. anserina seems to be a centralized network in its first phase of growth, the branching process dramatically and quickly increases the network complexity, meaning that the number of possible paths between the different regions of the thallus become very high. The centrality of the network is then lost, since the information can be exchanged without going through a center. Moreover, a partial destruction of the initial center (the ascospore location) does not affect the ongoing development of the network 1  www.nature.com/scientificreports/ three other nodes that can be distinguished by their chronological growth : there are indeed one earliest node, or mother node, and two daughter nodes. An illustration can be found in Fig. 2. Podospora anserina network can be distinguished from the four major types of networks previously mentioned. If this network was not a dynamic network as it was defined above, at a given point in time, P. anserina network could be locally related to a hierarchical and modular tree where modularity (existence of a basic architecture) is important, the heterogeneity (variance of the distribution of the number of links) is around zero and whose randomness (probability of creating links between the nodes) would be very low 15 . It is therefore possible to model the growth of P. anserina as different networks linked together by causal relationships and spatial constraints.
To implement this simulation we made the following assumptions: (1) the environment is supposed to be homogeneous (in terms of resources, light, and other exogenous constraints ...) and invariant over time; (2) the growth conditions are standard (temperature, pressure, etc.) and invariant over time; (3) there is no limits of deployment in space and time; (4) the initial state is defined as the germination of a single ascospore. These assumptions are not restrictive. Except for hypothesis (3) they correspond to the experimental conditions.
Basically, a simulation can be viewed as a two-step process. First is the generation step. It corresponds to the theoretical growth of the network (i.e. the theoretical law). This step is built from assumptions on probability laws  Table 2 for more details).

Figure 2.
Definitions used in this work. The thick line is the mother hypha, while the thin lines are the daughter hyphae. We consider two types of three degrees vertices. V 3ℓ are lateral branchings, connected with the angle θ ′ . V 3 are apical branchings on which two branches are connected: an operating and an exploratory branch with respective angles θ o and θ e . V 1 (resp. V 1ℓ ) are one degree vertices (i.e. apexes) coming from V 3 (resp. V 3ℓ ). www.nature.com/scientificreports/ and their parameters. Second is the detection step. It corresponds to the observation of this growth. This step is also built on probability laws but only relating to the observation process. It must reproduce as well as possible the entire experimental process. Details can be found in section of Observations of "Methods". This two-step process is justified by the fact that we never observe or measure what is "true" but only its convolution by the whole observation chain which can possibly generate biases, and necessarily has some resolution.
Finally this simulation allows us to propose a prediction on the distribution of apical branching angles.
Generation. The aforementioned considerations on the kind of network generated by P. anserina lead to choose a growth model in the form of a full binary tree (or proper binary tree) 16 . In this model, each branch (hypha) is driven by its tip (apex) that in turn is allowed to divide into a pair of sub-branches. This process is called apical branching in this paper, in line with the dynamic growth vocabulary. When apical branching occurs, the distribution of branch lengths and that of the angles between newly created branches with respect to the mother branch follow differentiated probability laws. An apical branching process is characterized by the emergence of a three bodies vertex, denoted V 3 , which is the point of intersection of the 3 branches and two one body vertex, denoted V 1 , the apexes of the branches (Fig. 2). Of course, for the duration of binary tree growth, both number of V 3 , N V 3 , and of V 1 , N V 1 , depend on the instant of observation and are related by where N ger is the initial number of branches emerging from the ascospore, or is large enough (or equivalently that the duration of growth is large enough). In the following, we will assume this condition is always verified. The length of the branches is a random variable assumed to be correctly described using a Gamma law which is a continuous and positively skewed distribution. Moreover, other laws of probability can be obtained by changing parameters of the Gamma law: with x the length, k the scale factor, θ the shape factor, and Ŵ(k) the Euler function. Using a particular choice of parameter values for k and θ , the Gamma law allows to fall back to (1) the exponential law ( k = 1 and θ = 1/ ), (2) the Maxwell-Boltzmann's law ( k = 3/2 and θ = 2a 2 ), (3) the χ 2 law ( k = ν/2 and θ = 2 with ν the number of degrees of freedom, d.o.f). Consequently, the Gamma law offers a lot of freedom to generate the length of the branches.
We considered a Gaussian probability law for the angles: N(x|x 0 , σ ) with x the angle, x 0 the mean and σ 2 the variance. Because it was empirically observed that apical branching leads to two different angles, the parameter relating to the mean, x 0 , is different for the two branches coming from a V 3 . The variances for the angles are chosen according to experimental data. The position of the branches in relation to the direction of the mother branch-to the right or to the left-is managed by a probability (i.e. Bernoulli law). So, it is impossible to force a particular direction for one or both forms in order to break the chiral invariance of the growth.
We call operating branch the branch coming from a vertex V 3 which has the most important angle (the angle is noted θ o ) with respect to the segment that connects the vertex to the vertex from which it came (mother vertex) and we call exploratory branch the branch with the smallest angle (the angle is noted θ e ) (see Fig. 2). This nomenclature is justified by the fact that the two branches do not have the same scope in terms of network growth: the branch showing the smallest angle tends to perpetuate the direction of the mother branch in order to explore the environment and then to capture new resources, while the other branch occupies space inside the network in order to optimally draw resources from the neighborhood.
Experimental observations show another type of branching process giving rise to another type of three bodies vertices: interior vertex leading to the creation of lateral branch. A branching occurring on a segment bounded by two pre-existing V 3 vertices, or a V 3 and a V 1 vertices can give a new lateral branch at a later moment during the growth process. To distinguish between these branches we call them lateral (operating) branch and the three bodies and one body vertices attached to them are denoted respectively V 3ℓ and V 1ℓ (see Fig. 2).
The creation of a lateral branch is also driven by a probability law. This law is unknown (as well as its parameters) leading to an important difficulty in implementing the phenomenon. However, we choose a power law based on two considerations: (1) the duration that separates the observation from the creation of the two vertices bounding the possible lateral branch (2) the length between the pre-existing vertices. In this case, the probability of generating a lateral branch as a function of the length and time parameters could be written: where x is the length of the branch (i.e. the length between the two vertices), x 0 a length parameter ( x 0 ≥ x ) (i.e. the greatest length generated by the probability law, typically the mean, kθ , plus three times the standard deviation, √ kθ ), α a parameter of the power law, typically α ∼ 2 and p 0 a scale parameter in the range [0, 1] to set the probability, typically p 0 ∼ 1/2 . The duration between the observation and the creation of the two vertices is not apparent in Eq. (2) because this probability is recomputed after each generation of newly generated branches (of hypothetical duration).
The location of the lateral vertices V 3ℓ between two vertices follows a uniform law. However we implement a censorship zone near the vertices so that the uniform distribution does not have the vertex-vertex distance as a support. Empirically, we see that the angular distribution between the segment connecting the two vertices and the lateral operating branch is about the same as the angular distribution of the operating branch. Therefore we choose the same law with the same parameters. Likewise, the length of the branches thus created follows the www.nature.com/scientificreports/ same probability law as that of the exploratory branches and operating branches (with different parameters). The position of the lateral branch in relation to the direction of the mother branch-to the right or to the left-seems to be in the direction of the curvature of the hyphae. Also, the chirality for these branches can be broken by a probability parameter. The control of these laws and their parameters is a challenge which is raised by the study of the calibration of simulation to the data (see the paragraph Calibrations... in "Methods"). The curvature between any two vertices is not taken into account in the present work for the following reasons: (1) the curvature can not be related to anything particular during the growth process (to a first approximation); (2) the chiral symmetry can be broken from the angular distributions if one wishes to obtain a global curvature effect. In this sense, the chirality breaking of the lateral branches takes into account the curvature of the segments between the points formed by the triplet, grandmother, mother and daughter-two adjacent segments; (3) from the experimental data it can be seen (not shown) that mean curvature (i.e. the curvature over the different fractions of the length between two vertices) is locally zero.
Because the experimental process introduces a time dependency through the reconstruction of network images and because we want to study the growth dynamics of the network, we define a time between vertex generations. At this step the time is given by number in generation agency between daughter-vertices and mother-vertex. For the lateral branches the time is given by the relative moment-in terms of generation-at the creation of the branch. However, this theoretical time will be scaled by the data (see calibration section) and we will take into account a possible growth interruption during apical branching process as well as differentiated apex velocities of exploratory, operating and lateral branches. For the simulation step, however, these velocities are given by a Gaussian probability law with average values and a small shape (standard deviation) and we assume these velocities are constant during growth.
To sum up, the generation time found is certainly not the growth time of the fungus. However, this time is necessarily proportional to the growth "true" time. The constant of proportionality must be obtained from calibration on the data.
Simulation can start with as many branches as desired. We make this number coincides with the empirical observation of the ascospore germination that shows three initial branches with an angular separation of approximately 2π/3 . We then apply the growth as a binary tree for each of three branches as described below. Figure 3 gives an example of the final state of growth.
In order to make the results of the numerical simulation comparable to the experimental observations, we introduced to the simulation output the same bias and uncertainties. Details of the procedure are given in the paragraph Observation... of "Methods" section. Parameters of the simulation are then tuned to reproduce the main experimental observations, for both time scale and space scale. In time, the typical growth dynamics of the apex number and the total length growth 2 a t is regained, including the respective rate exponents, as can be seen in Fig. 8. In space, the inertia tensor of the V 3ob vertex distribution were extracted. The eigenvalues of the diagonalization of this tensor allowed us to compare the global geometry of the numerical thallus to the experimental www.nature.com/scientificreports/ one, including sphericity (shape) and spreading (size). A detailed procedure of the complete calibration process is shown in the part Calibration... of "Methods" section.
Prediction of angle distribution. As noted above, when a V 1 vertex is converted into a V 3 vertex (an apical branching) the two "daughter" branches have very distinct angles. We can then make the hypothesis that the permanence of this behaviour is related to very deep reasons. There are undoubtedly biological reasons, but we can hypothesize that the source of this phenomenon is due to the fact that the fungus in a homogeneous environment has a clear advantage in occupying the largest surface in order to capture the maximum of resources from the environment and thus optimizes its growth. The surface occupied by the fungus can be written: with L the total length of the network, S bio the physical surface, or total area, of the fungus and S the surface of the intersections overlap (because the branches are flat tubes with constant width). For given growth duration and total length of the network, the surface occupied by the thallus is all the more important as the overlapping surface is reduced, i.e. the number of intersections is small. It is possible thanks to the simulation to test this hypothesis in order (1) to verify that indeed the largest occupied surface requires two quite distinct angles, and (2) to predict the optimal couple of angles. Figure 4 shows the logarithm of N V 3i (t, L, θ o , θ e ) �S(t, θ o , θ e ) as a function of θ o and θ e for a fixed time (i.e. generation index). Building upon the previous assumptions, we can make the following predictions, the wide angle is close to 80 • and the small angle is close to 15 • . The uncertainties on theses angles are estimated to be about 10 • because of the toy model used here (the fixed length of the branches at a fixed time about nine generations).

Experimental approach.
Podospora anserina is a coprophilous filamentous ascomycete, a large group of saprotrophic fungi that mostly grows on herbivorous animal dungs and plays an essential role within this complex biotope in decomposing and recycling nutrients from animal feces 17 . P. anserina has long been used as an efficient laboratory model to study various biological phenomena, especially because it rapidly grows at a rate of 7 mm/day on standard medium, it accomplishes its complete life cycle in only one week, leading to the production of ascospores, and it is easily usable in molecular genetics, cellular biology and cytology 18 .
Experimental set-up. An experimental device allowing to solve the dynamics of the local and global growth of the complete hyphal network of P. anserina directly on a Petri dish from an ascospore and over a period of ∼ 15 h in a controlled environment has been previously developed and described 12 . We did use of this setup to carry out three complete and independent series of images of the thallus growth, named (1), (2) and (3) thereafter, under standard growth conditions with M2 culture medium, see 19 for more details, at a temperature of 27 • C.
Direct measurement of angles. After the standard binarisation and vectorization process described in 12 we extracted angles θ o and θ e formed by an apical branching (see Figs. 2, 5) using direct reading on the network picture. At this stage, the thallus network is in its exponential growth phase, with approximately 1500 apexes and a length of 600 mm (see "Methods" for details, and more specifically Fig. 10; Table 2). We focused our analysis on www.nature.com/scientificreports/ apical branching. By observing the dynamic branching process real V 3 vertices can be segregated from the geometric vertices V 3i due to overlapping intersections. The comparison between images acquired before and after the branching allows the operator to extract a real V 3 vertices collection without any ambiguity. The procedure of angles extraction is shown in Fig. 5. A circle with a radius of 50 pixels (approximately five hyphal diameters) is centered on each V 3 vertex. The coordinates of the points of intersection between the circle and the three hyphae are manually recorded. For each apical branching, the extension of the mother hypha in the direction of growth determines two angles with both emergent hyphae, a small angle and a wide angle. From the complete network, 66 V 3 were extracted from each experiment. We then considered two populations, respectively composed of the collection of the small angles θ e and of the wide angles θ o (Fig. 5). The uncertainties on each angle value δθ = 4 • is estimated from the diameter of the hyphae. Figure 5 depicts the distributions of the measurements of θ e and θ o . Bin width reflects the uncertainties δθ . For the clarity of the representation, θ o defines the positive direction of rotation for each V 3 . For each extraction two populations are clearly visible. Assuming these populations follow Gaussian laws we extracted the respective estimators of the mean θ and width σ such as N(θ i |θ, σ ) using likelihood maximisation. The uncertainties on angle measurements is assumed to follow a centered Gaussian law with a standard deviation of δθ , such as N(θ i |0, δθ) . Building on the very tiny correlation between the estimators, we assumed cov(θ ,σ ) = 0 for the covariance matrices of the respective estimators. We finally compute an average estimate from this collection of i.i.d sets for each of the two angles, |θ o | = 70.9 ± 2.4 and σ o = 12.0 ± 2.4 , |θ e | = 9.3 ± 2.4 , σ e = 8.7 ± 2.4.
Since we must take into account the systematic uncertainties δθ = 4 • , the final results for the angles are: θ o = 71 ± 5 θ e = 10 ± 5.
GIS automatic method for automatic angle measurements. In addition to the work on direct measument angles, which focuses on a set of selected apical branching, we have implemented a Geographic Information System (GIS) automatic method for detecting, in a well-developed thallus, all angles for both apical and lateral branching. While being a non-selective method, it nevertheless offers a global information of branching process, with a high number of angle measurements. GIS geoprocessing was performed from the skeletonized thallus. Details of the skeletonization process can be found in "Methods" section.
The angles were calculated using the following procedure: discs (or buffers) of 5-pixel radius were centered on apexes ( V 1 ) and nodes V 3ob . All buffers are then splitted by the thallus skeleton. We tested many radius lengths in order to improve the method and we finally settled upon a smaller size radius than those defined for the direct measurement method (different by a factor of 10) due to the high density of line segments (i.e. the thallus) in certain portions of the image. Indeed, the latter would have generated additional splits of buffers by non directly connected line segments, which inevitably would have induced a bias in the measurements of disc areas and therefore of angles. The portions of discs are then converted into their equivalent angular counterparts. The results we obtained are not significantly different from one method to another. Details of this geoprocessing is shown in Fig. 6. A total of 1186 vertices have been detected, with 372 V 1ob and 814 V 3ob , respectively shown by red and black points in Fig. 6, leading to about 3000 angles.
Four populations of angles can be distinguished in Fig. 6. The peak observed at 0 • is an artifact due to a distance less than 10 pixels between two vertices. The peak at 360 • corresponds to V 1ob vertices that do not partition the circle. Finally peaks at roughly 50 • and 200 • correspond to the distribution of the three arcs of the circle partition. www.nature.com/scientificreports/ This pdf has been fitted using three Gaussians (see Fig. 6) of which one is only defined by a normalisation factor (the sum of the 3 angles must be 2π ). The χ 2 is very good for about 30 d.o.f. The estimated parameters are respectively 106 • ± 1 • and 137 • ± 3 • for the means and 8.4 • ± 1 • and 15 • ± 2 • for the widths of the interest gaussians (the third Gaussian is fully constrained because the sum of the angles must 2π).
There is only one geometric scenario to define θ e and θ o when we have 3 angles and when we know the hierarchy between these angles. The right scenario gives the angle values: θ e = 43 ± 6 • and θ o = 75 ± 5 • , taking into account the uncertainties on the angle constraint.

Discussion
In a previous paper 12 , we characterized the hyphal network expansion and densification of the filamentous fungus P. anserina-in particular through the evolution of the total length and the number of nodes and apexes-on a standard culture medium assumed to be homogeneous and optimal for the fungal growth. We have presented in details the experimental process implemented both from a biological and physical point of view. Here, we present a two-dimensional simulation of the fungal growth that allows us to better characterize some growth patterns of the fungal network. In particular, we predicted that thallus growth is driven by a specific angular branching process, optimized for the exploration and exploitation of the environment. In this work, a calibration of the simulation was first carried out from the collected data. To this end, we focused on the two main processes that determine the growth pattern of the fungal network-i.e. apical growth and branching-to propose a growth model in the form of a binary tree, meaning that each branch is divided into two sub-branches, as previously described 16 . Then, in the generation process, the main features of this model are that lengths and angles follow two different probability laws and that growth has to be considered dynamically. The fungal network growth has then been basically defined as an apical branching process characterized by intersection points between mother/ www.nature.com/scientificreports/ daughter branches ( V 3 ), and apexes of branches ( V 1 ). The growth model was incremented by lateral vertices ( V 1ℓ and V 3ℓ ) that can occur on a new daughter branch at a later moment of the growth. Such a process is in agreement with the usual description of fungal network formation 20 . Namely, two types of branching are allowed to occur according to their location on the hyphae: apical branching which is the emergence of two branches from a hyphal tip (named here V 3 ) and lateral branching from the sub-apical part of one hypha (named here V 3ℓ ). The proposed model was calibrated using the parameter values obtained from our experimental data 12 , namely the number of vertices for different times and the spatial geometry of the fungal network. The quality of experimental observations is limited by geometric intersections (overlapping). Obviously this limitation does not apply to the numerical simulation from which an estimation of the number of geometric intersections can be derived. Moreover, it unambiguously constructs the objective (and unobserved) history of the development of the thallus which is inscribed in the position of "real" V 3 . This simulation also allows an exhaustive study of chirality breaks during growth and/or at the origin of development. The most important feature of this model is the occurrence of two distinct angles, a wide angle and a small one at the daughter branches, when a V 1 vertex transforms into a V 3 vertex. This pattern of angular distributions at the apexes is not the one commonly used by the simulation attempts previously described. For example, Du et al. 7 has developed a lattice-based system for modelling mycelial growth of Postia placenta with an apical branching defined as the emergence of two branches that symmetrically develop with respect to the previous tip direction. Nevertheless, in our case this typical angular distribution allowed us to describe the apical branching process as the simultaneous emergence of two distinct types of branches. First is an exploratory branch, which slightly deviates (small angle) from the initial trajectory of the mother branch and could thus be involved in the centrifugal exploration of new territory, and then in the extension of the colony. Second is the operating branch, which shows a larger angle and deviates strongly from the initial trajectory of the mother branch. This type of branch could thus contribute to the exploitation of neighboring territories and then to the densification of the network. Such a pattern is in agreement with previous studies in which a clear distinction was made between leading hyphae at the edge of a fungal colony which grow into new territory in search for food, and the hyphae behind the colony edge that interconnect to form a three-dimensional network optimized to extract nutrients from the surrounding medium 4 . Then, leading hyphae, with an almost-linear trajectory from the center to the front of the thallus could correspond to the iterative development of exploratory branches from the apexes, contributing in turn to the extension of the thallus.
We hypothesized here that in a homogeneous environment, the fungus thallus growth is suited to occupy the largest surface in order to capture the maximum of resources and thus optimize its growth. Based on our hypotheses, our model predict that the largest occupied surface requires two quite distinct angles, a wide angle, close to 80 • and a small angle, close to 15 • (see Table 1). In previous studies, the measurement of branch angles has rarely been reported, and it remains in general limited in scope preventing the development of robust statistics. For example, Abd-Elsalam et al. 21 showed that morphological characteristics of Rizoctonia solani included a right-angled branching that could be used in identification of some isolates. Simonin et al. 22 carried out the measurement of angles on about a dozen leading hyphae of N. crassa and showed an angular distribution ranging between 50 • and 90 • . More recently, Du et al. 23 led a comparative plot of the angular distribution, both for apical and lateral branching on the thallus of P. placenta. They showed that this distribution is not affected by the age of region, or by the branch type and that branch angles remained approximately at 80 • , close to the right angle, which appeared to maximize the area covered by the colony.
From new sets of experiments, two approaches have been used in order to extract the branching angular distributions, (see Table 1). One is based on the extraction of the angles formed by an apical branching ( V 3 ) using a direct reading on the image, the other one derives from a GIS automatic method that allows for a global detection of angle measurements (apical and lateral branching). The main point of attention for the development of the first approach was the choice of the diameter of the circle allowing to measure the intersections between branches and thus to measure the branch angles themselves. A circle with a diameter of 50 pixels (approximately five hyphal diameters) was retained. This diameter was chosen in order to decrease the uncertainty due to hyphal width and to decrease the biases due to variation in the orientation of the apexes with time, independently from the initial branching. Branch angles extraction from approximately 200 apical branching allowed us to clearly identified the small and wide angle populations and to numerically estimate their corresponding values. The agreement is striking when compared to the optimal angular distributions that emerge from the prediction. However, the main disadvantage of this approach is that it requires the intervention of an operator who must select the apical vertices to be analyzed. This can be tedious, introducing a bias in the choice of apical branching, and necessarily result in a limited number of selected apexes. However, these results are confirmed by a second approach, complementary to the first one, that is global and based on GIS automatized detection of branch angle measurements. Here again, two populations of angles have been detected, with angular distributions that regain correctly the occurrence of a small and wide angles on apical branching. www.nature.com/scientificreports/ Comparison of numerical values leads to several methodological remarks concerning the wide and small angles. Wide angles were found in good agreement with direct measurement and prediction, within one standard deviation for the two angles. Small angles were found significantly higher for the small angle between GIS measurement and direct measurement (and simulation). Differences observed for the small angle measurements between the direct approach and the GIS approach could have several sources. Firstly, the choice of the vertices of the direct measurement could be slightly biased (the radius of the calculation of the angle is important which selects a certain type of vertices). Secondly, the radius of the circle delimiting the branching area for the GIS method, is smaller from the latter one by a factor of 10. This radius was chosen in order to optimize the statistics and avoid artefacts, notably in particularly dense areas of the thallus (high quantity of matter and high number of close nodes). However, this choice could introduce a bias on the result due to the linearization intrinsic to the method which is particularly sensitive for short lengths. Lastly, another bias could be due to the fact that no discrimination between lateral and apical branching could be established with the GIS based method, leading to take into account both cases. Moreover, there must be in the set of vertices allowing the reconstruction of the angles by the GIS method, a remainder of vertices V 3i . However, at three standard deviations, i.e. over 99.9% , the results remain consistent.
Overall, from both prediction and experiments we can conclude that, in our conditions, the process of angular branching allows P. anserina to occupy as much surface as possible, and then to explore and exploit its environment in the most optimized way. All these considerations must be replaced in an in vivo context, in which the fungus is usually in competition with other organisms to occupy the colonized area and the use of available resources. So, we will consider extending the calibrated model to simulate the mycelial growth in heterogeneous environments or under various constraints, as apical branch angles respect their intrinsic property but the extension of emerged branches is influenced by external factors. Likewise, since it is possible to calibrate the simulation in time and space on the experimental data, a study of the different speeds (biological, creation of V 1 , V 3 , V 1ℓ , V 3ℓ , group velocity, etc.) is conceivable and controllable. A more ambitious simulation is under development. It will allow the joint analysis of local curvature in relation to branch length or lateral branch emergence.
To conclude, our understanding of the coordinated growth and behavior of hyphae inside the fungal network is still in its infancy. In this context, our work contributes to show that a reasoned reductionism on the constraints of growth, both experimentally (two-dimensional growth) and mathematically (toy-model approach) is a powerful tool to describe relevant behaviour of the mycelial growth allowing for exploring multiple hypotheses and constraints rapidly.

Methods
Observation of the numerical thallus growth. The result of the measurements (or in an equivalent way, of the observations) always corresponds to a convolution of the reality, the probability density function (pdf) of "true value", with the acquisition chain (the pdf of resolutions). In order to make numerical results comparable to the experimental observations, it is then necessary to introduce the same bias and measurement uncertainties. In what follows, observations may refer to experimental or numerical data. There are different interferences which confuse the values of the quantities of interest that may depend on the entire detection chain which goes from the shooting to the post processing of the digitized image.
• First, the experimental observation process is not able to distinguish "true" V 3 vertices (i.e. from fission process or equivalently branching, V 3 and V 3l ) from geometric vertices (i.e. geometric intersection of two branches because of two-dimensional observation). The number of these geometric vertices depends on the time of observation (because of the creation of lateral operating branches and operating branches). Indeed, the probability of obtaining branch intersections increases with the density of branches which also increases with time. We distinguish two situations: (1) either the intersection is purely geometric (the experimentally observed vertex may result from the superposition of two branches. These vertices are noted V 3i because even this kind of vertex shows four branches. Usually, network reconstruction based on experimental images in this situation leads to two vertices with 3 branches in close proximity among the four branches constituting the intersection. In the second situation (2) it gives rise to a merger. We know that a process of hyphal fusion (i.e. anastomosis) can occur in P. anserina. These vertices are noted V 3o despite the fact it still is a vertex with four branches. The anastomosis phenomenon arises following an unknown probability. However, the number of vertices resulting from the hyphal fusion process is certainly very low compared to the number of V 3 + V 3l and V 3i vertices. Even if it is marginal, this hyphal fusion phenomenon allows a more efficient information transition within the network and it also makes it possible to be freed from a dependency in the centrality of the network. It is possible to qualitatively study, via simulation, the gain of a hyphal fusion phenomenon by building the adjacency matrix 13,14 and by iterating on it so as to connect all the vertices to each other in a minimum time (or distance). The presence of these vertices blurs the information on the "true" number of vertices V 3 , however, since N V 3 ≈ N V 1 and N V 3l = N V 1l for a long duration, we get a very reasonable order of magnitude of N V 3 + N V 3l . This limitation of experimental observation forces us to abandon the idea of working on historical markers of the development of a thallus. Indeed, if the vertices V 1 give us an image of the current state of development, the vertices V 3 -resulting from the V 1 which are only intermediate stagesconstitute the deep history of growth. In order to remedy to this corner case, we developed a method to clearly identify geometric intersections. • A second effect is the resolution. On the one hand, the optic we use obviously has a finite resolution. On the other hand, there are experimentally limits in the reconstruction of the images. We take these resolutions into account by convoluting the position of the different vertices positions and branches with a Gaussian www.nature.com/scientificreports/ resolution. Although taken into account in the simulation of the detection, this phenomenon is not of capital importance in the observation of the macroscopic characteristic quantities of the growth of the network. • The third effect is also due to the finite resolution of the optic. A vertex V 1 or V 1ℓ can be very close to a branch.
In this case it is impossible to distinguish it from the branch. This situation generates a vertex that behave as in the case of geometric intersections. We do not need to distinguish them from the vertex coming from intersections (i.e. V 3i ), and so we note them V 3i . It is easy to see thanks to the simulation that this phenomenon only becomes important when the density of hyphae is important, therefore for very long growth times. • Finally we also take into account the possibility of having ghost vertices on the images coming from the optics or from the reconstruction of the images. We note this vertices V 3g or V 1g . This type of vertex is minimicked by specific processing during image reconstruction. Also, its number should be marginal and does not depend on time. Here again, we can experimentally verify that this process is marginal.
The number of numerically observed vertices, N V 3ob and N V 1ob , is therefore the sum of the different contributions of the effects listed above, while dominant parasitic vertices are mostly of geometric nature. It is difficult to build a theoretical model which allows to obtain this number even if we know that this number grows very significantly with the time of observation (typically exponential growth). The simulation makes it possible to obtain this number.
Calibration of the numerical simulation. In the simulation, there are several parameters that need to be calibrated on the data. In order to do that we present here two sensitive variables.
Time scale. To scale the time, we use the number of vertices V 1ob and V 3ob . It is much more complicated to find the law that governs the number of V 3ob vertices than that of V 1ob because of the presence of V 3i whose law is unknown. On the other hand, because there are more intersections between branches than there are apexes near a branch, one would expect to obtain a larger contribution of branch-branch mergers, even if the merging process remains marginal (the V 1 contribution in N V 3o ). So, the calibration must be carried out on the vertex V 1ob (t) and we will check that the number of vertices V 3ob (t) follows the empirical law, see Eq. (7), from the data. The number of V 1ob (t) is formally: For reasonably long time and with a quality control procedure, we can neglect some terms of the sum, i.e. N V 1g (t) is very small compared to other numbers of vertices. The number of hyphal fusion apexes, ǫ 1 N V 3o (t) , is a second-order corrective term because it depends on anastomosis which is a marginal process. The number of apexes-branches intersections ǫ ′ 1 N V 3i (t) depends on the resolution of the optics and the reconstruction which can be controlled (it is always possible to degrade the resolution in order to check this number at a given time). However, we cannot effectively neglect this component even if the vertex density is not very high (we will see on the data observed experimentally that it occurs beyond a certain growth time). Also the calibration on the experimental data we use can only be carried out in a growth time domain.
The term N V 1ℓ (t) depends on the number of branches at the time t − δt (with δt the time between two observations in experimental data and index of generation in simulation) and the full probability of producing a lateral branch before this time. The number of branches at the time t − δt is . Finally because we are dealing with the growth of a binary tree, N V 1 (t) is know to be N V 1 (t) = 2 at with a > 0 , where N V 1 is a scale parameter for the time. Using these approximations the number of V 1ob (t) is then: with δt/t ≤ 1 and with p the mean fraction of lateral branches for a branch (which is a function of the probability to create lateral branches) which is constant by hypothesis. A is a constant amplitude which takes into account p and the fact that if δt is perfectly defined, the origin of the times is not: there is a "latency" time between the start of germination and the first observation. Indeed, we know that for t = 0 we must obtain 0 ≤ N V 1ob (0) ≤ 3 . We can then rewrite this number of vertices as N V 1ob (t) = A ′ 2 a(t+t 0 ) and in this case 0 ≤ A ′ ≤ 3 and t 0 is the "latency" time. So, A ′ 2 at 0 = A.
We fit experimental data with this law in order to get the 2 parameters (see experimental approach) then we fit N V 1ob (t) coming from the simulation with the same law. In order to calibrate the simulation time we get the relation (A2 at ) data = (A2 at ′ ) sim and compute t ′ with the hypothesis that t ′ is a linear function of t.
In Fig. 7 (left) , we plot the uncertainties for one standard deviation obtained from the fit parameters of data (grey area). We scale the time (and compute one standard deviation for this time from the uncertainties on the parameters of the 2 fits carried out) and plot the points coming from the simulation (the black cross, the uncertainties on the numbers of V 1 coming from simulation follow a Poisson distribution).
The law for the number of V 3ob is N V 3ob (t) ≈ 2 at n c n t n because N V 3 have the same law as N V 1 since N V 3i is the dominant term of N V 3ob and we assume a regular increase with time for this number. In practice, the series can be limited to the first 2 terms (check the first negligible order in the fit). Grey area in Fig. 7 (right) show one standard deviation giving from the fit parameters of data. If the scale factor obtained from V 1ob is correct, the points from the simulation for the three bodies vertices, must be inside a side band from the data.
The result shows a good agreement. We conclude that (1) we have correctly defined the time scale factor, and (2) binary tree growth is the right model for dynamic growth of P. anserina. However, it is certain that the law of growth of the fungus we use, A 2 a t ( a > 0 ), diverges when time tends towards infinity. This behaviour is certainly www.nature.com/scientificreports/ not permissible. Also, the law actually followed by this growth should rather be of the form N V 1ob (t) = 1 b+c 2 −a t with a, b, c real positive parameters. Using this law, we can see that if at is not too big and if c ≫ b , the law becomes N V 3o ≈ c2 a t which is the law of binary tree. So, the law we have considered in this analysis is valid for not too long (i.e. t − t 0 < 15 h) but not too short time (i.e. t − t 0 > 3h).
Space scale. With a growth model in the form of a binary tree, it is easy to obtain an observable which makes it possible to calibrate the distances in the simulation. For a given generation g, the number of vertices V 1 is where the factor three comes from the number of branches of the germination. The number of segments between two vertices ((V 3 , V 3 or ( V 3 , V 1 )) is N seg = 3 (2 g−1 − 1) , so that the average length of the network can be written �L� = �ℓ� N seg where ℓ is the average length between two vertices. The ratio, r 1 , of the average length of the network to the number of V 1 in the network is then r 1 = �ℓ� (2 − 2 2−g ) so that this is constant for g ≫ 1.
Because of lateral branches we need to correct this ratio. By construction, the number of lateral branches is proportional to the number of vertices segments. So, r 1 ≃ �ℓ� (2 − 2 2−g ) 1+p 1+2p for g ≫ 1 and where p is the mean fraction of lateral branches for a branch. r 1 is again a constant for g ≫ 1 . The index of generation g is linked to time that was calibrated previously.
So, for a time not to small we must show on data a constant for r data 1 . This is the case and we can estimate this constant on the data so as to calibrate the species on the simulation.
In Fig. 8 (left), we plot the uncertainties for one standard deviation giving from the fit parameters of data (grey area). We scale r sim 1 by the relation r sim 1 = r data 1 for a time not to small (using the scaled time obtain before) and plot the ratio coming from simulation (black cross). The uncertainties on the numbers of r sim 1 are obtain by considering a Poisson hypothesis and relative uncertainties for L tot equal 5%.
In Fig. 8 (right), we show the ratio r data 3 coming from the fit (grey area is one standard deviation) and we plot the r sim 3 for the simulation with the space scale factor and the time scale factor obtained by v 1ob . Note that it is impossible to define the law of r data/sim 3 because of the geometric intersections. However we know that this ratio must decrease as a function of time since the number of geometric intersections increases strongly and independently of the total length of the network. The agreement for r 3 between experimental data and simulation is convincing.
Another test is to check if the spatial geometry of the network has the same geometric behaviour between the experimental data and the simulation as a function of time. For this purpose, at a given time, we construct the inertia tensor of the vertex distribution for V 3ob (resp. V 1ob ): with x 0 , y 0 the mean position of the vertices cloud and x n , y n the position of the vertices.
Diagonalisation of this tensor gives two eigenvalues (and two eigenvectors, which are the main axes of the vertices cloud). These eigenvalues are ordered as 1 ≥ 2 and we define an indicator of the geometry of the vertex cloud: s = 2 2 1 + 2 (called sphericity). If 1 ≫ 2 then s → 0 , which means that the vertices cloud has the shape of a elipse, if the eigenvalues are degenerated, s = 1 , which means that the vertices cloud has the shape of a disk. www.nature.com/scientificreports/ In the data the shape of the vertex cloud is concerted over time ( s(t) = constant ). Its value depends on the initial conditions, i.e. the way the spore was generated. So the algebraic value of the constant is not relevant and only the shape of s(t) is relevant.
We check if the simulation gives a constant sphericity with scaled time and space. Figure 9 gives s(t) for the simulation and for the data of the V 1ob distribution in the plan. The values for the uncertainties are obtained by a bootstrap method 24 . The sphericity for the V 3ob vertices are also constant for data and simulation but because of the V 3i it is difficult to give a right interpretation.
We check whether the direction of the main axes of the vertices cloud move as a function of time. In the case of experimental data this direction is constant even s ∼ 1 . The same applies in the case of the simulation.
Observation of the thallus growth. After the standard binarisation and vectorization process described in 12 we extracted and adjusted the following quantities in order to calibrate the simulation process. Let N V 1ob be the number of observed V 1ob , N V 3ob be the number of observed V 3ob and L the total length of the network. We assume the uncertainties associated to the vertices count to be Poisson. The total length L is about N V 1ob + N V 3ob segments of average length ℓ . An estimation of the uncertainty associated with the measure of the total length L was derived as σ L = N V 1ob + N V 3ob �ℓ� . We assume the uncertainty associated to the acquisition time t is  www.nature.com/scientificreports/ half the sampling period δt . The growth of such a microorganism follows basically a five steps timeline: (1) a lag phase, (2) an exponential phase, (3) a deceleration phase, (4) a stationary phase, and (5) a decline phase. We study only the tree first steps. We previously showed that the temporal growth of N V 1ob , L and N V 3ob , to a first approximation, can be modelled by an exponential growth. This of course can only be valid locally (i.e. for short growth times) and from a purely descriptive perspective. Indeed, to be consistent with the binary tree model we made use of a base-2 exponential function to describe the growth of these quantities. Thus the following expression may be used: where X i stands respectively for N V 1ob , N V 3ob or L and t is the time. τ is the characteristic growth time, t 0 is the temporal offset corresponding to the transition between the end of the lag phase and the first observation. For the experimental data, there is necessarily a time offset between the first recorded image and the ascospore germination. In order to compare the experiments with one another, this time offset should be taken on a per experiment basis because it varies according to each experiment. X 0 i are the respective parameters N 0 and L 0 , corresponding to the growth onset given the temporal offset. N 0 V 1ob , N 0 V 3ob and L 0 values are respectively expected to be approximately 3, 1, and in the range 10-20 hyphal diameters. Equation (7) shows a non-linear behaviour of the fit parameters. Consequently we made use of the following procedure: We excluded spurious data, i.e. for time greater than 15 h and for time smaller to 2h. The value of t 0 was manually adjusted in order to obtain simultaneously N 0 V 1ob , N 0 V 3ob and L 0 in the range of the respective expected values. We assume that the uncertainty on t 0 corresponds to 2δt . Note that t 0 gives access to the lag phase duration. With t 0 as a fixed parameter the least squares fit is linear. The grayed areas in Fig. 10 show the results of the fit (for one standard deviation) on the considered data range. The points are the data with their uncertainties. Solid black lines represent the best fit parameters (see Table 2). The respective grey shadowing wraps the data range used to perform to statistical fit, while its thickness quantifies the associated uncertainties to one standard deviation. www.nature.com/scientificreports/ The growth rate exponents τ and the measurement of specific quantitative parameters for the three experiments are summarized in Table 2 and Fig. 10, respectively.
Because these 3 measurements are i.i.d. we assume the final results for the pertinent parameters of the growth are those shown in Table 3. Note the difference between the doubling time τ V 1ob and τ V 3ob . This is a straightforward consequence of the contribution of the geometric vertices in the three body vertices collection. Conversely, doubling time for the one body vertices τ V 1ob and total length τ V L are found to be in agreement.

Skeletonization stage and vertex detection for GIS automatic method. Several methods for
generating centerlines, or skeletons, from polygon features have been described in the literature. We can cite Voronoi diagram based method, medial axis transformation algorithm or Delaunay triangulation method 25,26 . Among those, an automatic Geographic Information System algorithm, ESRI Polygon to Centerline-available online thanks to the ESRI 27 plateform-appears to be efficient for modelling skeletons of elongated polygonal features 25 . The Polygon to Centerline GIS tool used in this work, based on the creation of Thiessen polygons, is described as follows: (i) The first step is to create a large number of vertices from a thallus polygon feature (methodological aspects were discussed in 12 ) through an adaptive densification method for improving the centerline quality. Vertices are then converted into point features (red point features in Fig. 11a). Table 2. Summary of the growth rate exponents for N V 1ob , N V 3ob and the total length L extracted from data shown in Fig. 10. The χ 2 values were all found to fall into the range 1-5.  Table 3. Final results of the experimental growth rate exponents, from data shown in Table 2. www.nature.com/scientificreports/ (ii) Then the Thiessen polygons may be created, each of them contains a single point input feature (i.e. a vertex). We can observe a high number of Thiessen polygons located both inside and outside the fungal thallus, and of course only Thiessen polygons located inside the thallus area interest us. Thus, we extract polygons that overlay the thallus by using the "clip" geoprocessing tool (represented in beige in Fig. 11b), which are then converted into linear features. All linear segments that are connected to the thallus single polygon are removed to constitute the centerlines, represented in purple in Fig. 11b. (iii) At this step of the treatments, we can see that a high number of short purple linear segments have to be removed (or trimmed). We used the ESRI "Trim Line" tool, which consists in removing in our case all purple lines that do not touch another line at both end points. The result of this treatment is shown in Fig. 11c (the orange plus green segments). (iv) A small number of branching artifacts still need to be removed (the orange ones). The hyphal width being relatively homogeneous, orange centerlines are removed by calculating a morpho-mathematic erosion (creation of a three pixels 'negative' buffer from the thallus envelope): all the segments that touch the three pixels wide surrounding area (grey shaded area in Fig. 11c) are trimmed and only the green segments are kept to form the skeleton. It should be specified that this method is efficient to produce centerline features and for working on V 3ob vertices, but it slightly alters the position of the apexes ( V 1ob ). (v) We then used successively two GIS tools to dissolve all centerline segments into an unique linear feature and to generate a feature class of points for locating the vertices (nodes and apexes). At this processing step, each V 3ob vertex is materialized by three superimposed point features when each V 1ob vertex corresponds to only one point feature; this makes it possible to distinguish V 3ob vertices from V 1ob vertices. Finally, V 3ob features are dissolved to produce one feature per vertex (instead of three) (see Fig. 6.