Calculating Hyphal Surface Area in Models of Fungal Networks

—Filamentous fungi grow efﬁcient nutrient transportation networks which are highly resilient to attacks by grazers. Understanding them may beneﬁt the design of human-built networks, where such properties are sought after. We recently developed a mathematical model that improved previous 2-dimensional studies by representing the space in a 3-dimensional face-centered cubic lattice. While the model focused on structural aspects (hyphal orientation, branching, and fusion), these are closely tied to functional aspects, that is, the handling of nutrients. In this paper, we reﬁne our previous model by modelling the hyphal network as a set of cylindrical tubes connecting spherical junction points, and calculating the exact local hyphal surface area. In further development of the model, this will allow the reﬁnement and incorporation of existing nutrient consumption models—in particular, how nutrients are used for turgor maintenance at a particular network location.


I. INTRODUCTION
Fungi, like all kingdoms of life, are economically and ecologically critical. Fungi have two modes of growth, and a number of mathematical models have been developed to better understand and control their growth. In yeasts, a minority of fungi, cells separate after they divide. In filamentous fungi cells do not physically separate and are kept together in filaments called hyphae. These hyphae extend at their tips and branch, forming a network, the mycelium, that penetrates the environment and absorbs resources. In this paper, we focus on the mathematical modelling of filamentous fungi.
A large body of research has been devoted to modelling filamentous fungi, as covered by Prosser in the late 70s [10], summarized by Kotov and Reshetnikov in 1990 [9], or reviewed more recently by Davidson in 2007 [3]. Models can be intuitively divided based on the scale they focus on, since fungi can cover hectares due to indeterminate growth 1 while their basic machinery operates on the micrometer scale [3]. In this paper, we focus on the micrometer scale, which is of particular interest when attempting to understand how the function and structure of fungi give rise to specific patterns of growth. Models at this scale can be further divided into three types [1]: macroscopic models, focusing on quantities at the colony-level (e.g., overall growth rate, biomass density), microscopic models, explicitly representing every hypha, and intermediate models which represent quantities over several hyphae and are classically reaction-diffusion models. In order to precisely understand the complex coupling of structure and function in fungi, the most accurate mathematical abstraction is offered by microscopic models.
L Bakker, A Poelstra, Calculating Hyphal Surface Area in Models of Fungal Networks A. Contribution of the paper Boswell and colleagues developed a model at the micrometer scale in which hyphae are explicitly represented [2]. This model was two-dimensional, and thus created only planar fungal networks (i.e. there cannot be a hypha on 'top' of another). We recently introduced a three-dimensional model to better reproduce the structure of fungi [4].
The move to a three-dimensional model makes it particularly challenging to compute the surface area of the mycelium, and the cost of turgor maintenance. Fungi work to generate osmotic forces, consuming resources to generate turgor pressure, and drive water uptake and the bulk flow of cytoplasm towards the growing hyphal tips [7], [8].
Hyphae regulate local osmotic pressure (and thereby water uptake) by transporting ions from the environment into the hypha or vice versa. The cost of ion transportation will be a function of the osmotic pressure gradient and the surface area exposed to this gradient. Thus the local turgor maintenance cost depends on the available hyphal surface area.
This cost needs to be taken into account to determine which hyphae have a surplus of nutrients available for growth or branching. The key contribution of this paper lies in calculating the exact local hyphal surface area in the model, which improves on both [2] and [4] and can be further used to increase the accuracy of fungal network models.

B. Organization of the paper
In Section II, we explain how the physical space is discretized in our model. As our focus is on the cost of turgor maintenance, we refer the interested reader to [2], [4] for the rules detailing the orientation, branching, and fusion of hyphae in the model.
Having described a discretization of space, in Section III we fix the remaining free parameters (e.g., lattice width) and choose our coordinate systems. We see that the lattice width is restricted by the cross-section radius of the hyphae. Our hyphae are constructed by connecting cylinders and spheres. In Section IV we compute the surface area of these hyphae by determining how much of the surface of the spheres is not covered by cylinders, and how much of the surface of the cylinders is not covered by other cylinders.

II. DISCRETIZED SPACE
Dividing the space into chunks involves two questions: dimensionality and degrees of freedom. Firstly, the dimensionality asks whether the fungus should be modelled as spreading in a slice of soil (2 dimensions), or spreading as it naturally does in a volume (3 dimensions).
Limiting a model to two dimensions means that hyphae cannot be overlapping, thereby forcing the structure to be artificially planar. Therefore, we investigate the 3dimensional case, which allows for overlapping hyphae. Secondly, the degrees of freedom must be sufficient so that an extending hypha can change direction with a realistic angle. As illustrated in Figure 1(a), if the cells were square then the hyphae can extend to the cell in front (angle of 0 • ), the cell on the right or on the left (angle of 90 • ). Thus, the hyphae would either keep its direction or change it by 90 • . An angle of 60 • was suggested as a more accurate description [2], leading to the hexagonal cells depicted in Figure 1(b). In fact, 60 • is the smallest allowable angle in a uniform discretization of 2D (or 3D) space.
Obtaining an angle of 60 • between neighbours is straightforward in two dimensions: the hexagonal grid shown in Figure 1(b) is the only space-filling configuration of equal-sized cells that achieves this. In three dimensions, there are several arrangements that produce 60 • angles between neighbours. Such arrangements fall under two schemes: Hexagonal Close Packing (HCP) and Face Centred Cubic (FCC) arrangement [5], [6]. Of these, FCC is the only viable candidate (Figure 2(a)), since the asymmetry in HCP can prevent a hypha from growing straight ahead, which is its most common direction ( Figure 2(b)). FCC is an extension from the 2D hexagonal grid, in which horizontal hexagonal layers are stacked on top of each other but displaced slightly to allow for a tight packing. This extension can be intuitively understood by induction. . This is not the case for the local neighbourhood in a Hexagonal Close Packing (HCP) arrangement, which is symmetric only along some edges between the central cell and a neighbour (b). As a consequence, a fungal tip that grew onto the central cell from the right below (red line) would not be able to grow exactly straight ahead.

III. DISCRETIZATION IN 3D
Assume there is a discrete model simulating the growth of a mycelium in three dimensions. The model must take into account the amount of resource needed for hyphal wall maintenance for each discrete point in the lattice to determine where excess resource permits growth. In the FCC lattice chosen in Section II, each vertex of the lattice is assigned the part of the hyphal surface that lies within its Voronoi cell. Therefore, knowing the cost of hyphal wall maintenance requires computing the surface area represented by each vertex of the FCC lattice. Formally, given a subgraph G of an FCC lattice with a sphere of radius r at each vertex and a cylinder of radius r and length ∆x at each edge, we want to compute the external surface area of the union of all these shapes.

A. Notation
From here on, we will use spherical (θ ∈ (−π, π], φ ∈ [− π /2, π /2]) coordinates, where θ represents the angle between a vector and the positive x axis and φ represents the angle between a vector and the positive side of the x, y plane; and cylindrical coordinates (θ ∈ (−π, π], z ∈ R), where θ is the same as before and z represents the height above the x, y plane.

B. Adequately discretizing the space
In order to 'charge' maintenance costs to individual vertices, the hyphal surface area needs to be computed for the Voronoi cells associated with these vertices. This puts restrictions on the size of the spheres and cylinders that are put at vertices and edges, respectively. If their radius r is too large in relation to the Voronoi cell size ∆x, a fraction of a hypha's surface area may be assigned to a cell that it does not logically belong to (see Figure 4). Therefore, ∆x should be chosen such that a hypha between two vertices, say v i and v j , falls entirely within the boundary between the associated pair of Voronoi cells and does not 'invade' other cells.
Theorem 1. Hypha connecting two vertices v i and v j falls entirely within the associated Voronoi cells when Proof: The intersection of a hypha (represented by a cylinder) with the boundary face between v i and v j is a circle of diameter 2r. The centre of this circle is at the intersection of edge e ij and the boundary face. Without loss of generality, pick any diameter of the circle and put two vertices, v k and v , as close as possible to either end of the diameter (see Figure 5), but at distance ≥ ∆x from The cylinder between vi and vj needs to fit through the boundary between their Voronoi cells, which runs between the centroids of vivjv k and vivjv . v i and v j . This maximally constrains the boundary size, which is the distance between the centroids of v j v i v k and v i v j v .

C. From two to three dimensions
In 2 dimensions the surface area (i.e., surface length since it is 1-dimensional) within a particular hexagonal Voronoi cell v i with a set H i of edges occupied by hyphae is equal to ∆x + πr when |H i | = 1, or where H i is a clockwise ordered circular list, i.e.
. This is illustrated in Figure 6. In three dimensions, the first term of (1) needs to be adjusted only to account for the extra dimension: 2 · |H i |·πr ∆x 2 . The second and third term generalise to three dimensions, although a one-pass scan (as in (1)) does not suffice anymore. The method we will use for both the portion of the sphere to be added and the portion of the overlapping cylinders to be subtracted are similar. For a generic sphere (cylinder), construct the refinement of all overlaps with (other) cylinders. Then for every one of the 2 12 possible neighbourhoods, mark the appropriate areas in the refinement as 'exposed'. Store the resulting areas in a lookup table with the neighbourhood (represented as a binary string) as lookup key. The second and third terms of (1) are to be replaced by lookups in these tables. The construction of the refinements used to create these lookup tables is detailed in the next section.

IV. AREA CALCULATION
In Section III, we described a discrete model simulating the growth of a mycelium in three dimensions. In order to compute the surface area contained in a Voronoi cell in the lattice, we required 'maps' of which parts of the surface of a central sphere is covered by which neighbours' cylinder, and of which part of the surface of a neighbour's cylinder is covered by other neighbours' cylinders. In this section, we describe these 'maps' and compute the required areas.

A. Central Sphere
For every hypha that is connected to a vertex v i , a hemisphere of that vertex's sphere is covered, and therefore removed from the surface area. What remains Vertex labelling congruent with directions ϑ k ; opposite directions (−ϑ k ) omitted.  is a hemisphere, lune, spherical triangle or spherical quadrangle (e.g., Figure 7- Figure 10). Each neighbour v j divides v i 's sphere into two hemispheres along a circular boundary that is perpendicular to e ij . The refinement of all divisions associated with a neighbour of v i is constructed as follows. The six directions ϑ k along which neighbours are placed are numbered in Figure 7. We denote the direction opposite of ϑ k (mirrored in c) as −ϑ k .
The intersection of the boundaries corresponding to neighbours v j and v j are given by ±(ϑ j × ϑ j ). 1 . These intersections lie along one of seven directions ϕ k : The actual refinement is constructed by connecting the intersection points ϕ k along the associated circular boundaries, to produce a polygonal partition of the sphere's surface area.
Theorem 2. The partition of the surface area of v i 's sphere as induced by its neighbours' cylinders consists of 24 triangles of equal size.
Proof: We first observe that the partition will consist of polygonal regions, where each polygon vertex is an intersection point ϕ k .
At ±ϕ 3 , ±ϕ 6 and ±ϕ 7 , we have the intersection of two circles, producing four π/2 angles each, 24 such angles total. (These angles can be calculated explicitly or derived from symmetry of the lattice.) Each pair of 1 Alternatively, the intersections could also be derived from the convex hull of the neighbour positions (each ϕ k corresponds to the centre of one of its faces), or from the Voronoi Diagramme of the neighbour directions on the sphere (each ϕ k corresponds to one of the Voronoi Diagramme's vertices). The refinement is the Delaunay Triangulation of the ϕ k . L Bakker, A Poelstra, Calculating Hyphal Surface Area in Models of Fungal Networks these intersection points is separated by a circle, so that each polygon in the partition has at most one of these intersections as a vertex (and therefore at most one interior angle of extent π/2). Therefore the partition has at least 24 polygons.
In total, the intersection points are responsible for 72 angles, which are the interior angles of the polygons in the partition. This gives an upper bound of 24 polygons, which is achieved when every one is a triangle. This, combined with our earlier observation that we have at least 24 polygons, tells us that the partition consists of exactly 24 triangles.
Finally, since the area of a spherical triangle is entirely determined by its angles, which are π/2, π/3, π/3 for each triangle, every triangle has the same area.
To determine the topology, we notice that each cylinder directly touches six intersection points. For any of the remaining eight points, say, φ, either (a) every triangle with φ as a vertex is covered, or (b) none of them are. Which case applies for a given vertex is geometrically obvious. Using this, we determine the topology of the refinement, and create Figure 12. This figure also encodes the relationships between a neighbour and the areas of the sphere it covers: each neighbour lies at the center of a circle that its cylinder would cover, e.g. ϑ 2 lies at the centre of the circle through ϕ 2 , ϕ 7 , −ϕ 5 , −ϑ 5 , −ϕ 2 , −ϕ 7 , ϕ 5 and ϑ 5 . A cylinder in direction ϑ 2 would cover all areas inside this circle, and a cylinder in the opposite direction (−ϑ 2 ) would cover all areas outside this circle.

B. Hyphal Cylinder
The cylinder that is put at every hypha e ij not only covers part of the surface of the spheres belonging to the cells it connects, but it also covers (and is covered by) parts of the surface of other cylinders connected to these cells. The cylinder associated with another hypha e ik , j = k connected to one end of the cylinder intersects it along an elliptical boundary ζ k that is defined by the plane through e ij × e ik , e ik × e ij and e ij + e ik . As was done for the central sphere, the refinement of all divisions associated with an edge e ik needs to be constructed.
Cutting the cylinder associated with ϑ 1 along the plane θ = 0, we see that these intersections between the cylinder associated with e ij and those associated with e ik , j = k correspond to sine waves ( Figure 11).
Specifically, if a hypha ϑ i has coordinates (θ, φ), the intersection of its cylinder with that of ϑ 1 will have the form which follows immediately from converting to ϑ 1 's cylindrical coordinates.
The refinement of these boundaries is a subdivision of the cylinder into 30 triangles and 6 quadrangles, of 11 different sizes. These sizes can be calculated as sums of definite integrals over the sine waves associated with the edges. The symmetrical nature of the lattice supports that areas which appear equal actually are equal. The only feature of this diagram not geometrically clear is that the 3-way intersections marked • are true; i.e., they are not just three 2-way intersections very close together.
It suffices to verify one 3-way intersection. This is because the diagram will look identical no matter which hypha we use as a base, up to relabelling-by replacing ϑ 1 by ϑ 3 , ϑ 4 or ϑ 6 , respectively, the other three intersections marked with • can be positioned at the •.
We see that the leftmost • is a true 3-way intersection, since when x = 0, the three curves have the value ϑ 3 : tan ( π /6) sin ( π /2) = 1 /3 ϑ 4 : tan ( π /4) sin cos −1 2 /3 = 1 /3 ϑ 6 : tan ( π /3) sin cos −1 8 /9 = 1 /3 All intersection points lie along one of 26 directions. Below are the coordinate specifications for the intersection points in the interval   Fig. 11. The cylinder associated with ϑ1, when cut in two along the vertical plane θ = 0, illustrates that the area calculation consists entirely of intersections of areas below a sine wave. The 11 areas numbered in red are the only ones for which the area needs to be calculated. The others follow by symmetry. Like in Figure 12, intersection points i j ϕ are labelled. To distinguish the labels from those used in Figure 12, the intersection is labelled using the indices of two of the arcs that cross (ζi, ζj). i j ϕ : θ , z −5 −6 ϕ ∼ −ϕ 5 : The actual refinement is constructed by connecting the intersection points along the associated circular boundaries. By a similar process to the one used to create Figure 12, we can determine the topology of the refinement, and create Figure 13. This figure also encodes the relationships between a neighbour e ik and the areas of e ij 's cylinder it covers: each neighbour lies at the center of a circle that it covers, e.g. ϑ 2 lies at the centre of the circle through ϑ 5 , ϕ 1 , ϑ 5 , ϕ 1 , ϕ 16 , ϕ 20 and ϕ 22 . A cylinder in direction ϑ 2 would cover all areas inside this circle. Unlike with the sphere, a cylinder in the opposite direction (ϑ 2 ) would not cover all areas outside this circle; rather, ϑ 2 has its own circle (through ϑ 3 , ϕ 1 , ϕ 13 , ϕ 18 , ϑ 3 , ϕ 17 , ϕ 11 and ϕ 1 ). −ϑ k or −ϕ k denote ϑ k or ϕ k reflected in the origin, respectively. It has been projected from ϕ1 onto the plane θ = π /2, and rotated counterclockwise by π /2 to emphasise the symmetry.
V. CONCLUSION We have computed the surface area for paths in a FCC lattice, which model hyphae growing in a mycelium. This improves on earlier models in that it is 3D and can therefore model non-planar mycelia, the typical case in nature. By using a hexagonal lattice rather than a cubic one, we allow our hyphae to grow in straight lines or to change direction by as little as 60 • , giving a realistic range of motion while still being computationally feasible.  Fig. 13. Topology of the refinement of cylinder boundary intersections (ϕ k ) induced by the 11 neighbours (ϑ k in blue) of the central cell. −ϑ k or ϕ k denote ϑ k or ϕ k reflected in the origin, respectively. It has been projected onto the plane from above ϑ1.