 Research
 Open Access
 Published:
A framework for an organellebased mathematical modeling of hyphae
Fungal Biology and Biotechnology volume 2, Article number: 5 (2015)
Abstract
Background
Although highly desirable, a mechanistic explanation for the outstanding protein secretion capabilities of fungi such as Aspergilli is missing. As a result, a rational and predictive design of strains as cell factories for protein production is still out of reach. The analysis of the secretion apparatus is not only hampered by open issues concerning molecular cell biological processes, but as well by their spatial fragmentation and highly dynamic features. Whereas the former issues are addressed by many groups, an account of the space and timedependent processes, which is best done by means of mathematical models, is lacking. Up to now, mathematical models for hyphal organisms mainly focus on one of two extremes. Either macroscopic morphology, such as pellet or mycelium growth, is addressed, or a microscopic picture is drawn predicting, for instance, the form of a hyphal tip. How intrahyphal transport and organelle distribution works, however, has not been tackled so far mathematically.
Results
The main result of this contribution is a generic modeling framework to describe the space and timedependent evolution of intracellular substances and organelles. It takes intrahyphal, passive and active transport of substances into account and explains exponential and then linear length growth by tugordriven uptake of water. Experimentally observed increasing concentration levels of organelles towards the tip can be well explained within the framework without resorting to complex biological regulations. It is shown that the accumulation can be partly explained by geometrical constraints, besides a necessary deceleration of the active transport velocity. The model is formulated such that more intricate intracellular processes can be included.
Conclusions
Results from steadystate experiments are easy to be interpreted. In a hyphal network, however, new branches are produced at an exponential rate. Moreover, passive and active transport processes give rise to a spatial distribution of organelles and other cytoplasmatic constituents inside hyphae. As a result, most of the data obtained in experiments will be from a nonsteady and space dependent state. A quantitative and mechanistic explanation of the processes occurring will only be possible if these dependencies are taking into account while evaluating experimental findings.
Background
The ecological and technical relevance of fungi is outstanding. They are integrated in most ecosystems, act as detrimental agents for plants and humans, decompose waste materials, and are exploited in the synthesis of valuable products [1, 2], to name just a few areas in which they play a major role. Their most striking feature is polarized growth and branching which leads to more or less dense mycelia or pellets [3, 4]. Concomitantly, Aspergilli such as A. niger, A. oryzae and A. terreus have astounding capabilities to secrete interesting enzymes, mainly through the apical region [5]. A rational design to obtain modified strains as optimized cell factories, however, is still limited by the incomplete picture of their growth, production and secretion machinery. By the very nature of living cells all occurring processes are highly dynamic and the behavior of a cell does not only depend on the actual stimuli but what had happened to the cell in the past. For fungal organisms, interpretation of physiological data is even more challenging. Besides the compartmentalization of biological functions in distinct organelles, space and timedependent distributions occur. This relates to organelles, other cytoplasmic compounds, and stimuli in and around a mycelium, and, therefore, impede the deduction of meaningful knowledge and hypotheses. All of which could be addressed in the context of mathematical models.
So far most mathematical models for hyphal organisms respecting morphological features focused on macroscopic processes in dense mycelia and pellets. Here, the interplay between nutrient transport by diffusion and spacedependent growth was addressed, see for example [1, 6–18] and references therein. A second group of models tries to predict the geometric appearance of rather small mycelia [19–25]. In these models, very little detailed biological information is needed, if used at all, to give rather realistic pictures. The results of simulations shown in Figure 1, as an example, just use the information that the apical growth rate in three dimensions is constant and that septa and branches are formed when a critical length of a hyphal compartment is obtained. Including a random growth direction in the simulations gives a visual impression which is hard to distinguish from real photographs. As for the first group of models, these models will not help in deciphering the secretion apparatus. This is true as well for socalled morphologically structured models, which do not even account for the spacedependency, see, e.g., [7, 21, 26].
A more detailed account of microscopic features of individual hyphae is given by a last group of models to describe, for instance, the shape of a tip, or the growth in length. A wellknown example for the first class of problems, which will be used in what follows, is given by BartnickiGarcia et al. in [27]. Here, the geometrical form of an apex is predicted with a simple model. It is based on a set of hypotheses how vesicles are transported ballistically from the Spitzenkörper to the wall. The model has been refined in future works to better account for the threedimensional shape of a tip, the way vesicles are transported to the wall by diffusion, by representing the cell wall as a flexible membrane, or by a better account of vesicle fusion with the cell wall [28–32].
Mathematical descriptions of the growth in length, as the second class of microscopic features, were given recently with two different approaches. In [33], the longrange transport of material in hyphae is depicted by a particle transport along a single, hypothetical microtubule extending over the whole length of a hypha. The amount of material reaching the tip of the hypha determines length growth. A changing velocity, however, is neither considered nor a movement of the microtubules with the cytoplasmatic flow. In contrast, [34] explain length growth of Phanerochaete velutina mathematically by a turgor driven intrahyphal flow towards the tip. In all these approaches a constant length growth rate is considered which is not true for the germ tube. Moreover, new branches of a mycelium very often show a lower initial velocity as well. As a result, and as a mycelium grows exponentially by an exponential production of new branches, a significant part of a mycelium will not be in a kind of quasisteady state which is assumed above. In major parts of a mycelium, organelles and intracellular substances have not yet reached their quasisteady state distribution, which might be important for a quantitative prediction of the growth of the mycelium. Likewise, if septa are closed and opened by Woronin bodies, intrahyphal flow has to stop or will resume resulting in even more complex situations.
In our former works [11, 19], we explained the initially observed exponential and then linear growth with the limitation by a hypothetical intracellular compound. We had to resort to a hypothetical compound at that time as neither for fungi nor for actinomycetes details about the mechanism were known. Especially for fungi, this situation has changed drastically in recent years. Molecular methods, bioinformatics and image analysis have provided us with a whelm of information if not give rise to a ‘Big Data Tsunami’ [35]. More specifically, for the processes addressed here, which are responsible for length growth and (product) secretion, much more is known today. Excellent recent reviews about growth in length of fungal hyphae are given, e.g., by [36, 37], and about secretion in [32]. The importance of turgor driven length extension is stretched by [38, 39] in a series of papers.
We therefore think that the time is ripe to try to condense at least a small part of the available knowledge in a mathematical model. This can form the basis to discuss hypotheses and to account for the effects of space and timedependencies in the interpretation of experimental data. The model structure derived in this contribution therefore serves two main purposes.

First of all, it represents a basic model structure, with which the initial exponential and then linear growth of a hypha can be described with a minimal amount of assumptions. This turgor driven evolution of the intrahyphal flow forms the ‘backbone’ for all other processes occurring in a hyphae, and, therefore, has to be considered first. To be more specific, besides the postulation of some kinetics, no further biological regulations will be introduced to describe the experimentally observed growth evolution. If this is possible, already simple physical transport processes combined with implicitly formulated regulations through kinetic expressions can be used to explain the observed behavior without resorting to complex biological mechanisms. This, of course, does not rule out such regulations which additionally may occur. If experimental evidence is given, these processes can be included readily.

Secondly, the model structure derived serves as a basis for future work when experimental data is interpreted and condensed in a mathematical framework. As an example, the distribution of vesicles in a hypha will be considered here which shows a distinct profile along the length of a hypha. Again simple physical arguments, mainly with respect to the active transport velocity and the geometry of the tip, will be enough to explain experimental data where a significant increase in concentration is observed toward the tip.
In the long run, such kind of models might help in answering questions raised in the endeavor toward a rational strain design. Examples are [5]: How many vesicles carrying proteins of interest can be used without interfering with vesicles for growth? Where are the bottlenecks in vesiclemediated protein secretion? How many proteins can be channeled through the secretory pathway in order to provide each protein sufficient time to become correctly folded? Extending this list of questions will naturally occur when a model is at hand.
The rest of the paper is organized as follows. After a problem statement in the next section, the model of BartnickiGarcia et al. is revisited to determine the volume and surface area in the tip region. This will then be used to correct experimental data. The general model is formulated next. As a first application, length growth by turgor driven water uptake is described. Extending the model with vesicles allows for a comparison with experimental data in the last section before the paper finishes with some conclusions.
Problem statement
What will be described below, will neither encompass the process of sporulation nor branching. An attempt to model germination can be found, e.g., in [40]. It is assumed that a very short part of a tip already exists, presumably from a mother hypha. In the simulations shown, depending on the boundary conditions, the mother hypha will either have no influence on the developing branch or it will supply material for growth. After some time, say \(t=t_t\), the retrograd end of the new hypha will posses a final radius of \(r_t(t_t)=R\) and the length of the tip will be \(L_t(t_t)=L_{tmax}\). For \(t>t_t\), the hypha will consist of this constant tip part with fixed length \(L_{tmax}\) and a seemingly growing distal region of length \(L_d(t)\), i.e., the overall length of the hypha will then be given by
Consequently, it will be assumed that, after tip completion, its geometrical form and size stay constant. Both assumptions together simplify the mathematical treatment significantly as apparent growth is associated to a subapical region where the cell wall becomes rigidified, probably by the action of crosslinking enzymes [37]. With this assumption, the complex processes involved in plasma membrane and cell wall synthesis actually occurring in the tip region are approximated. A more detailed model is given by [32]. A constant form and size of the tip region, on the other side, necessitates a constant radius of the distal part. This is in contrast to the model introduced in [34] where the authors correlate high internal volume flows Q with large diameters of hyphae. This will not be considered here for simplicity although the general model could account for it.
For a simpler numerical implementation, the tip and the distal part will be modeled separately with appropriate conditions to account for the connection between them, see Figure 2. In what follows, the variables related to the tip will be denoted by an index t whereas those related to the distal part with a d.
The general model will be obtained by a formulation of balance equations applied to an infinitesimal intrahyphal balancing volume of length dx, see the space between shaded areas in Figure 2. Terms of production, consumption and transport via the cell membrane, the consequences of an intrahyphal flow, and, finally, due to active translocation will be considered. Spacedependent uptake of nutrients could be included readily, but is not done here. Most importantly, a constant physiological and functional state is assumed along a hypha. If this is not the case, the model developed in this contribution has to be combined with approaches proposed, e.g., by Nielsen and Villadsen [18]. As the transport from and to the environment will be proportional to the local surface area and production and consumption rates will be given based on the local volume, these quantities have be determined first.
Geometrical model of the tip
For simplicity, the 2Dmodel proposed in [27] describing the form of a hyphal tip is used to derive expressions for the local surface area \(A(x_t)\) and volume \(V(x_t)\). Although the initial conjecture that the actual 3Dform of a hypha can be produced by a rotation of the solution of the 2Dmodel was corrected in [28], see as well [29], the simpler approach is used here. This is motivated by the fact that the actual differences in the forms obtained are small while the calculation of the 3Dform is rather involved.
In [27] it is proposed that: (1) the cell surface expands from materials discharged by walldestined vesicles, (2) vesicles are released from a postulated vesicle supply center (VSC), (3) vesicles move from the VSC to the surface in any random direction. Based on these propositions, they derive the following model.
If the VSC is located in the origin of an \((x_h, y_h)\)coordinate system, as used by BartnickiGarcia et al. [27], where the cartesian coordinate \(y_h\) equals the radius \(r_t\) in Figure 2, i.e., \(y_h=r_t\), the 2Dgeometry reads
In this coordinate system, the hypha extends to negative \(x_h\)values. The foremost point of the tip is at \(x_h=D\). With \(x_t=x_h+L_tD\), see Figure 2, the geometry of the tip is given by
in the coordinates used in this work.
For the calculation of the volume \(V(x_t)\) and surface area \(A(x_t)\), \(r_t(x_t)\) is locally approximated by a straight line of length \(l_t\) connecting \(r_t(x_t)\) and \(r_t(x_t+dx)\). Rotating this line defines the area and the enclosed volume.
Neglecting higher order terms, the infinitesimal surface area of such a truncated cone is given by
and the infinitesimal balancing volume is approximated by
The derivative needed in the former expression can be obtained from Eq. 3 by implicit differentiation, see Appendix A.
For a stationary observer, see Figure 3, the balance volume increases when the tip grows out of a considered section \([x_t, x_t+dx]\).
As it is assumed that the form of the tip stays constant, every point of the surface moves with the extension rate \(\dot{L}\) of the hypha, hence,
A negative sign is included as for a stationary observer the volume increases while \(d r_t/d x_t\) is negative.
The tip and the distal part must be connected without a step in radius, i.e., \(r_t(x_t=0) = R\). As a result, D cannot be chosen independently when \(L_{tmax}\) is fixed. From a biological point of view, fixing D might be a better alternative. This, however, would lead to a very long tip region which is ruled out here for numerical reasons. From Eq. 3, when \(L_{tmax}\) is prescribed, an implicit expression is given for D
For small \(L_{tmax}\), this value of D is an approximation of the real distance of the VSC from the apex. Figure 4 gives an impression of the obtained form of the tip, where some distal part is shown as well. The calculations were done with \(R=3.5\) μm, \(D=1.2\) μm, and \(L_{tmax}=15.6\) μm.
Experimental data
Although the main goal of this contribution lays in the derivation of a generic model structure, some comparisons with experimental data will be done.
Unpublished experimental data with vesicle measurements from Aspergillus niger is kindly supplied by F. Spanhoff, A. Ram and V. Meyer. They visualized the secretory vesicle concentration of individual hyphae by the intensity of the fluorescent RSNARE protein SynA using a Zeiss confocal microscope. Pictures were taken at an equidistant spacing of 0.2 μm along the length of a hypha for typically 10–12 layers (zstacks) across the diameter. Hyphae were obtained from the periphery of a fungal colony. Uncalibrated fluorescence data is obtained by adding up all intensity values for each zstack to obtain I. The mean value for 7 different hyphae is shown in Figure 5 as a function of the distance to the apex. The origin of the xcoordinate system used here is chosen \(D=1.2\) μm behind the apex, which coincides with the VSC used later. As the scanned volume in the apical dome is smaller than in subapical regions, the data is corrected here for this geometrical effect. Again the model of [27] is used to calculate a local volume, and, hence, from the intensity data, a local normalized measured concentration \(V_m\). As absolute information about the number of vesicles is missing, the data is additionally normalized by an arbitrary scaling factor of 35 to compare against simulation results later. The scaling factor is chosen such that in the subapical region a normalized concentration of approximately 1 is obtained. The corrected data is shown in Figure 6. Observe the higher ratio of maximal to subapical values in the corrected data compared to the given intensity values in Figure 5.
In the experiments described above, hyphal length was not measured for A. niger as a function of time for newly developing branches. Typically, an initial exponential growth will be observed followed by a linear one. As the model will be able to describe this, another set of data is used here for comparison. Experimental results from [41] are exploited. This rather old set of data was already used by us in [11]. In this former work, a much simpler model was proposed to describe length evolution. Using the data again, both approaches can be compared. Fiddy and Trince [41] measured the evolution of a primary branch of Geotrichum candidum extending out of an intercalary compartment just behind a septum. They observed a correlation of the decreasing length extension rate of this branch with septation occurring after some time in this branch. However, the extension rate of the primary branch continued to increase, despite septation, until a length of about 700 μm was reached. From Figure 3b in [41], a maximal extension rate of 2.5 μm/min can be estimated. The data will be given later together with the results of a simulation.
Generic model
To derive a generic model, a substance \(\mathcal{S}_i\), \(i=d,t\) is balanced in a segment extending from \(x_i\) to \(x_i+dx\), see Figure 2. In what follows, \(\mathcal{S}_t\) represents vesicles (\(\mathcal{V}_t\)), Spitzenkörper (\(\mathcal{K}_t)\), osmolytes (\(\mathcal{O}_t\)), etc. in the tip, i.e., \(\mathcal{S}_t = \{\mathcal{V}_t, \mathcal{K}_t, \mathcal{O}_t, \ldots \}\). Accordingly, \(\mathcal{S}_d = \{\mathcal{V}_d, \mathcal{K}_d, \mathcal{O}_d, \ldots \}\) will denote variables in the distal part. A radial distribution inside the hypha and diffusion in all directions are neglected. Diffusion in the xdirection could be included readily without complicating the numerical solution much. It would, however, make less sense for organelles. See [39] for a discussion of diffusion coefficients of different cellular components compared against the intrahyphal flow velocity.
The volumetric mass concentration of a substance \(\mathcal{S}_i\) is represented here by the very same symbol in the equations, i.e., \(\mathcal{S}_i\) is used to denote the concentration of the generic compound. Hence,
Balancing the mass \(m_{\mathcal{S}_i}\) of substance \(\mathcal{S}_i\) in a segment of volume \(V(x_i)\) with surface area \(A(x_i)\) and infinitesimal length dx leads to
with production (\({}_{\ldots p}\)) and consumption (\({}_{\ldots c}\)) rates \(\mu _{\mathcal{S}_ip}\) and \(\mu _{\mathcal{S}_ic}\), respectively. Transport (\({}_{\ldots t}\)) from or to the surroundings is modeled by \(\mu _{\mathcal{S}_it}\). Whereas the former reaction rates are given as a temporal change of mass per volume \(V_i=V(x_i)\), the latter is based on the local external surface area \(A_i=A(x_i)\) of the segment.
The last two terms in Eq. 9 represent intrahyphal flow in and out of the balance volume, i.e., flow through the shaded areas in Figure 2. As due to turgor pressure hyphae take up water from the surroundings, and as only the apical region can extend in the real hypha, an intrahyphal flow is set up. Hence, the volumetric flow rate \(Q_i(x_i,t)\) is both a function of space \(x_i\) and time \(t\).
Intrahyphal flow of a substance can be the result of the flow of the cytoplasm (cyt) transporting \(\mathcal{S}_i\). For other substances, longdistance transport is realized via an active (act) dislocation along microtubules. In the latter case, microtubules can be transported as well with the flowing cytoplasma, see [42, 43], resulting in a superposition of flow and dislocation velocities for \(\mathcal{S}_i\),
Alternatively, it can be assumed that microtubules stay fixed with respect to the cell wall and that all of \(\mathcal{S}_i\) is attached to them. Then, \(Q_{i,cyt}=0\) for this specific substance.
The splitting up of \(Q_i\) will now be used for a volume balance which is only affected by the cytoplasmic flow. As the cytoplasm is incompressible,
where \(\mu _{V_ip}\) represents the volume production, e.g., through turgordriven uptake of water from the environment through the local surface \(A_i=A(x_i)\).
To finally set up the generic model structure based on the balances given above, several steps are necessary which are detailed in the Appendix B:

1.
The last term of Eq. 9 and the second term of the righthand side of Eq. 11 are expanded in a Taylor series, neglecting all terms in \((dx)^n\), \(n \ge 2\).

2.
All equations are combined.
This leads to
with
Initial and boundary conditions, and \(Q_i\) will be specified below.
Generic model of the constant, distal part
For simplicity, it is assumed that the active translocation velocity along microtubules in the distal part, if it exists at all, is constant, i.e., for its gradient
Only at the tip, a deceleration will be considered later. Furthermore, after setting \(i=d\), and \(A_d= 2 \pi R dx\), \(V_d= \pi R^2 dx\), Eq. 12 reads
with
The fourth term on the right hand side can be interpreted as a dilution term due to intrahyphal flow. From Eq. 31, see Appendix, with \(\dot{V}_d=0\), an expression can be given describing the spatial evolution of the intrahyphal flow
As long as the volume production \(\mu _{V_dp} \ne 0\), Eq. 16 gives rise to a monotone increase of \(Q_d\) with \(x_d\), i.e., more and more fluid will be transported toward the tip.
The boundary conditions, \(Q_d(0,t) = Q_{dx0}(t)\) and \(\mathcal{S}_d(0,t) = \mathcal{S}_{dx0}(t)\), describe the information coming from a spore or from a branching site of a mother hypha. Whereas a mathematically ‘convenient’ boundary condition \(Q_{dx0}(t)=\) const. would make sense, as it describes an active transport, which can be zero as well, \(\mathcal{S}_d(0,t) = \mathcal{S}_{dx0}(t)=\) const. would be more difficult to justify biologically. This would mean that the mother hypha or spore would not change its value of \(\mathcal{S}\) irrespective of what is going on in the new hypha. This could only be explained by a source of \(\mathcal{S}\) of infinite strength. If the spore or the mother hyphae are not determined by a separate model,
is a better choice, at least from a numerical point of view. An initial condition for \(\mathcal{S}_d\) must not be specified as the simulation will start without a distal region.
Before specifying the individual production and consumption rates this generic model equation will be adapted to the nonconstantarea and nonconstantvolume case seen in the tip.
Generic model of the tip
Due to the nonconstant surface area \(A_t\) and volume \(V_t\) of the balancing volume the expressions get more involved. With Eqs. 4, 6 and 31 the resulting cytoplasm flow reads
The first term on the right hand side increases the flow due to turgor driven volume production \(\mu _{V_tp}\). The second term, however, as \(dr_t/dx_t\) is negative, decreases the flow to account for the volume needed for length growth in the tip region.
Using Eq. 12,
with
results. A plot of these \(\rho _t\)terms is given in Figure 7. In the distal part, the corresponding terms are constant, see Eq. 15.
Besides an initial condition \(\mathcal{S}_t(x_t,0) = \mathcal{S}_{tt0}(x_t)\), a boundary condition has to be specified. If a distal part is not yet formed,
is applied. If, on the other side, a distal compartment already exists,
Model of length growth
As the radius of the distal part and the form of the tip stay constant, \(\dot{L}\) is given by the overall volume produced divided by the area of the growing end. If the tip is still developing, this area is \(\pi r_t^2\). Hence,
When the tip is finished, with \(L_t=L_{tmax}\), it is assumed that length growth is realized by extension of the right hand side of the distal part with area \(\pi R^2\), see Figure 2. For a tip with constant geometry, this translates into a growth in length of the hypha of
Here, a hypothetical flow rate \(Q_{tp}(L_t,t)\) at the tip \(x_t=L_t\) has to be used. It is determined by the gross increase of the volume flow rate in the tip region
with a proper boundary condition for \(x_t=0\). Actually, the real flow rate in the tip region has to decrease towards the tip as it is ‘consumed’ everywhere in the tip to fill up the new volume formed due to tip growth. In the model, however, length increase is attributed to an increase in the distal part or left hand end of the tip before it is fully developed. Therefore, as volume is taken up in the tip region as well, the gross increase of the volume flow rate has to be known to determine \(\dot{L}(t)\). To distinguish this hypothetical from the real flow rate, it is denoted by \(Q_{tp},\) instead of \(Q_{t,cyt}\). The boundary condition at the left hand side of the tip, see Figure 2, is \(Q_{tp}(0,t) = Q_d(L_d,t)\) or 0 when only the tip exists.
Complete generic model
In summary, after specifying the individual production, consumption and transport rates, \(\mu _j\), and the change in active transport in the tip, \(\partial Q_{t,act} / \partial x_t\), the following equations have to be solved in the generic model assuming that the tip has already reached its maximal extension:

1.
Integration of Eq. 16 determines the intrahyphal flow rate at \(x_d=L_d\) which sets the boundary condition for Eq. 23.

2.
Integration of latter equation leads to the hypothetical flow rate at the tip, \(Q_{tp}(L_t,t)\),

3.
and, with Eq. 22, to the actual extension rate \(\dot{L}\).

4.
In a moving boundary framework, as L(t) grows, Eqs. 14 and 18 are solved to determine \(\mathcal{S}_i\), \(i=d,t\).
Initially, only the tip region exists. Hence, Eqs. 16 and 18 are not needed.
Modeling pressure regulation via osmolytes
We consider osmolytes which are responsible for maintaining a certain pressure and pressure gradient inside hyphae, see [38]. It is assumed that osmolytes are produced until a certain pressure is obtained for which intracellular sensors must exist. For the MAPK pathway, OS1 is discussed as a sensor in [39]. For simplicity, an intracellular substance called osmolyte \(\mathcal{O}_i\), \(i=d, t\), is introduced, which represents both the osmolyte, and, indirectly, the pressure. To obtain a mass flow toward the apex, its concentration must be higher in subapical parts. Using the equations derived above, \(\mathcal{S}_i\) is now replaced by \(\mathcal{O}_i\).
To start with the most simple model, it is assumed that osmolytes are not transported actively, \(Q_{i,act}=0\), and are not consumed or degraded, i.e.,
Furthermore, they are not transported over the cell wall, hence,
For the production of \(\mathcal{O}_i\), as a first approach, a logistic lawlike expression is used
with a maximal osmolyte concentration \(\mathcal{O}_{max}\). It has to be pointed out that these assumptions can be changed easily.
In the long run, for all osmolytes produced, water has to be taken up. If this process is fast, it can be assumed that, in a kind of quasisteadystate point of view, water is taken up proportional to the synthesis rate of the osmolytes. As the first one is formulated based on the surface and the latter based on the volume,
leads to an expression for the volume production rate \(\mu _{V_ip}\) in part i of the hypha
With Eq. 24, for the distal part,
and for the tip
follows. Different constants, \(k_2\) and \(k_3\), are introduced to possibly account for the fact that most water is taken up in the tip region, as hypothesized by [39]. This difference in water uptake velocity might be the result of the plasticity of the wall in the tip region. As the maximal extend of the tip, \(L_{tmax}\), considered here is fixed arbitrarily, this has to be observed during parameter identification and analysis of the simulation results.
With these kinetics and the generic equations derived in the last section, the models of the distal part and the tip can be formulated. They are omitted here for brevity. Furthermore, to reduce the number of kinetic parameters and to ease parameter identification, a normalization is done with \(o_i = \mathcal{O}_i / \mathcal{O}_{max}\) and \(q_i = Q_{i,cyt} / Q_{max}\), and \(o_i, q_i \in [0,1]\) for \(i=d,t\), \(\dot{L}_{max}=Q_{max}/(\pi R^2) = Q_{max} \rho _{d2}\), \(\theta _1 = k_1\), \(\theta _2 = k_2 \mathcal{O}_{max}\), \(\theta _3 = k_3 \mathcal{O}_{max}\).
For the numerical solution of the partial differential equations, the spatial coordinate is discretized equidistantly with a step size of \(\Delta x=0.2\) μm. The method of lines is applied for the equation describing the evolution of the osmolytes approximating the spatial derivatives by a firstorder backward difference operator. In the beginning, when only the tip exists, the left most discretization segment of the tip, see Figure 2, is allowed to grow in length according to \(\dot{L}\) until it exceeds a length of 0.3 μm. After the distal part is formed, its right most discretization segment takes over this task and grows accordingly until it exceeds a length of 0.3 μm. Then, this segment is split up into a segment of constant length (0.2 μm) and a growing one with an initial length <0.2 μm and the calculations are continued as before. The normalized flow rate \(q_i\) is obtained accordingly from Eqs. 16 and 23 exploiting the trapezoidal rule. Parameters \(\theta _i\) have to be chosen such that \(q_i \le 1\) is guaranteed.
As this work concentrates on the formulation of the generic model and not on a parameter fit or selection of appropriate kinetic expressions to describe a very specific problem, a simple approach was chosen to find kinetic parameters for simulation studies. Measurements performed by Spanhoff et al. were done with hyphae growing with approximately \({\dot{L}}_{max} = 3 \mu\)m/min. The parameters \(\theta _{1,2,3}\) are determined such that L(t) shows an initial exponential increase followed, after a transition, by a phase of constant growth velocity of approximately \({\dot{L}}_{max} = 3 \mu\)m/min. To this end, an optimization problem was formulated. In lack of real data for this first study, a ’desired’ evolution \(L_{des}(t)\) was determined to allow for an adaption of the \(\theta _i\)’s.
Simulation results are shown in Figures 8 and 9. Parameters used are \(\theta _1=0.3079\), \(\theta _2=0.2619\), and \(\theta _3=0.5032\) with vanishing flow from the mother hypha or spore. The initial tip starts with four discretization segments, i.e., a length of 0.8 μm, with exponentially decaying values of the initial, normalized osmolyte concentration \(o_{t,k}= [0.750\), 0.7351, 0.7206, 0.7063].
In the simulations given in Figure 8, growth is followed up to a length of about 170 μm, i.e., in the end more than 850 discretization segments are used. Due to production and consumption of the osmolyte, and due to length growth, the osmolyte profiles change dynamically over time. These osmolyte profiles, on the other hand, determine the overall volume production, see Figure 9, and, hence, length increase. As a result, all processes are highly interwoven.
When the rate of osmolyte production is increased by choosing higher values of \(k_1=\theta _1\), slower growth results, see left part of Figure 10. Here, as in the right part of this figure, the other parameters are as in Figure 8. Lowering \(\theta _2=k_2 \mathcal{O}_{max}\) has the similar, though less pronounced effect. An explanation of this feature can be seen in Figure 11 where the evolutions of the osmolyte and the normalized flow rate are shown for an increased value of \(\theta _1=0.3079 + 0.05\) compared to the case of Figure 8. A larger osmolyte production rate leads to higher values of o(x, t), and, hence, to a lower volume production rate \(k_i (\mathcal{O}_{max}  \mathcal{O})\), \(i=2,3\).
Using the same volume production parameters \(\theta _2=\theta _3=0.2616\) in the distal and tip region results in a much too slow growth corroborating different uptake rates. These simple studies show that an intricate relation exists between parameters and the processes described by the model. Predicting the outcome of parameter or kinetic variations is difficult. Therefore, numerical studies are indispensable.
In the first simulation study shown above, which will be used again in the next section, linear growth already starts at a length of approximately \(L=30\) μm. By another choice of parameters, the exponential phase can be elongated. In Figure 12, as an example, parameters are chosen such that experimental data obtained by [41] can be described by the model. In this case, an elongated exponential phase can be observed. Truly linear growth does not start before 700 μm as pointed out by the authors. As the primary branch considered in [41] originates from a parent compartment, the normalized flow rate coming from this compartment was fitted as a parameter as well. For simplicity, a constant value \(q_i(0,t)\) was assumed. A more detailed study should be performed, though, based on the general model introduced above, to better account for a variable influence of the mother compartment. Likewise, septation could be included. This would lead to a leveling out of the profiles inside the mother compartment and a variable supply of the new branch. Such a detailed study, however, is outside the scope of this contribution and has to be postponed to future work.
Modeling of vesicle distribution
As a first example of an organelle, vesicles are considered which have been already used above to get an idea of the form of the tip. They are, beside other functions, responsible for material transport to the growing tip. According to the model of [27], vesicles are used up in the apical part of a hypha. Therefore, the consumption rate in the tip will be \(\mu _{\mathcal{V}_tc} \ne 0\), while \(\mu _{\mathcal{V}_dc} = 0\) is assumed in the distal region.
As vesicles are produced inside hyphae, the transfer rate to the surroundings is zero in both parts, i.e., \(\mu _{\mathcal{V}_it}\), \(i=d,t\). For vesicle production \(\mu _{\mathcal{V}_ip}\), again a logisticlike expression is chosen. However, as vesicles are accumulated in the tip region, the expression is modified describing a production only as along as \(\mathcal{V}_i < \mathcal{V}_{max}\) and no production otherwise
Active transport is modeled by microtubules. It is assumed that all vesicles are bound immediately to microtubules. Therefore, the vesicle concentration can be calculated as the product of the local microtubules concentration, \(m_{MT}/V\), times the loading of the microtubules, \(m_{\mathcal{V}_i}/m_{MT}\)
Due to the discrete nature of microtubules, this is a rough approximation. Still, if the local microtubules concentration stays constant up to the apex, the active volume flow is simply given by
where \(u(x_i,t)\) and \(r_i\) represent the local transport velocity and hyphal radius, respectively. With this expression, the volume flow at the apex is zero which does make sense as \(r_t(L_{tmax})=0\). The change in flow rate over the length of a hypha is given by
In the distal part, with \(r_i=R =\) const.,
For a constant velocity u,
as assumed above. For the tip,
When the second term in the last equation is neglected, an increase in vesicle concentration, as seen in Figure 6, can only be observed directly behind the tip due to the decelerated transport coming from the geometrical influence (\(dr_t/dx_t < 0\)). This does not coincide with experimental findings. For this reason, as a first trial, it is assumed that u stays constant in most of the distal part at \(u=u_{max}\). From \(u_{max}\) it decreases linearly to zero in the foremost 15 μm. For the following simulation, \(u_{max}=35.62\) μm/min is chosen. Typical literature values for kinesin1 are in the range of 24–54 μm/min. For kinesin 7 in A. nidulans, [44] report 10 μm/min compared to the conventional kinesin of this strain with 120 μm/min.
The equations are normalized again with \(v_i = \mathcal{V}_i /\mathcal{V}_{max}\) and \(\theta _4 = k_4\).
Consumption of vesicles is assumed to be proportional to the actual length growth rate for which it is used, and proportional to the local vesicle concentration \(v_t\)
As the volume corrected measurements by Spanhoff et al. indicate a falling vesicle concentration only for the last three segments with a width of \(\Delta x = 0.2\) μm, for the rest of the hypha \(\mu _{\mathcal{V}_tc}=0\) is used.
A comparison of simulated data against the measurements introduced above is given in Figure 13. Additional parameters used are \(\theta _4=0.0407\), and \(\theta _5=2.2701\). The remaining parameters \(\theta _i\), \(i=1,2,3\), are as in Figure 8.
With these parameters the experimentally observed accumulation profile of vesicles can be reproduced. A keyenabler is the model assumption that the transport velocity must decrease towards the tip. This could be tested in future experiments.
Conclusions
As constantly new branches and septa are produced in a mycelium and a flow of cytoplasma towards the different apices occurs, quasisteady state concentration profiles of intracellular substances and organelles hardly establish. Moreover, due to these processes and due to active and cytoplasmatic flows intracellular components are not only a function of time but as well a function of space inside a hypha. A comprehensive and quantitative interpretation of experimental data of individual hyphae will therefore only be possible when these space and timedependent processes are taken into account. To this end, a generic mathematical model is proposed here which first of all describes turgor driven length extension. By this, an initially accelerated and then linear growth can be predicted as seen in microscopic experiments. A much simpler model, see [11], could be used to predict the length evolution though. That model, however, cannot to be extended so easily to describe other constituents of a hypha as it is done here. The turgor driven length extensions forms the ‘backbone’ of a generic model to study, for instance, the timedependent distribution of organelles and other compounds. These may be transported actively or passively towards the tip. Diffusion is not considered yet, but can be included readily. For a complete specification of the model, kinetic expressions have to be stated. In this contribution, very little effort is invested to study the effect of different kinetics and parameters, e.g., with respect to osmolyte or vesicle production and consumption. The emphasis is rather on the formulation of a generic model. Effects of different kinetics will be tackled in future works when more experimental data is available. However, even with simple kinetics chosen here it can be shown, for example, that the experimentally observed accumulation of vesicles near the tip can be explained. A crucial assumption to be able to do this is the postulation of a decreasing active transport velocity in the tip region. Without this, the experimentally observed apical accumulation of vesicles cannot be described in the chosen setting. The model structure can be readily extended to study the effect of different organelles and cytoplasmatic compounds. Before doing so, however, the numerical solution of the partial differential equations with a moving boundary should be revisited to hopefully decrease the computational burden. This was not done yet, as septation and branching have to be included in future works.
References
 1.
Papagianni M (2004) Fungal morphology and metabolite production in submerged mycelial processes. Biotechnol Adv 22(3):189–259
 2.
ElEnshasy HA (2007) Filamentous fungal cultures—process characteristics, products, and applications. Bioprocess Value Added Prod Renew Resour New Technol Appl 9:225–261
 3.
Metz B (1976) From pulp to pellet. PhD thesis, Technische Universität Delft, Delft
 4.
Kossen NWF (2000) The morphology of filamentous fungi. Adv Biochem Eng Biotechnol 70:1–33
 5.
Meyer V, Fiedler M, Mitsche BH, King R (2014) The cell factory Aspergillus enters the Big Data era: opportunities and challenges for optimizing product formation. Adv Biochem Eng Biotechnol. doi:10.1007/10_2014_297
 6.
Cohen D (1967) Computer simulation of biological pattern generation processes. Nature 216:246–248
 7.
Megee RD, Kinoshita S, Fredrickson AG, Tsuchiya HM (1970) Differentiation and product formation in molds. Biotechnol Bioeng 12(5):771–801
 8.
Edelstein L (1982) The propagation of fungal colonies: a model for tissue growth. J Theor Biol 98:679–701
 9.
Prosser JI (1995) Mathematical modelling of fungal growth. In: The Growing Fungus. Springer, pp 319–335
 10.
Tough AJ, Prosser JI (1996) Experimental verification of a mathematical model for pelleted growth of Streptomyces coelicolor a3 (2) in submerged batch culture. Microbiology 142(3):639–648
 11.
King R (1998) Mathematical modelling of the morphology of streptomyces species. Adv Biochem Eng Biotechnol 60:95–119
 12.
Lejeune R, Baron GV (1997) Simulation of growth of a filamentous fungus in 3 dimensions. Biotechnol Bioeng 53(2):139–150
 13.
Davidson FA (2007) Mathematical modelling of mycelia: a question of scale. Fungal Biol Rev 21(1):30–41
 14.
Krull R, Cordes C, Horn H, Kampen I, Kwade A, Neu TR et al (2010) Morphology of filamentous fungi: Linking cellular biology to process engineering using Aspergillus niger. Adv Biochem Eng Biotechnol 121:1–21
 15.
Celler K, Picioreanu C, van Loosdrecht MCM, van Wezel GP (2012) Structured morphological modeling as a framework for rational strain design of Streptomyces species. Anton Leeuw Int J G 102(3):409–423
 16.
Davidson FA, Boswell GP, Fischer MWF, Heaton L, Hofstadler D, Roper M (2011) Mathematical modelling of fungal growth and function. IMA Fungus 2(1):33
 17.
Hopkins SM (2012) A hybrid mathematical model of fungal mycelia: tropisms, polarised growth and application to colony competition. PhD thesis, University of Glamorgan
 18.
Nielsen J, Villadsen J (1994) Morphologically Structured Models. In: Bioreaction Engineering Principles. Springer, pp 229–269
 19.
Yang H, King R, Reichl U, Gilles ED (1992) Mathematical model for apical growth, septation, and branching of mycelial microorganisms. Biotechnol Bioeng 39:49–58
 20.
Meškauskas A, Fricker MD, Moore D (2004) Simulating colonial growth of fungi with the neighboursensing model of hyphal growth. Mycol Res 108(11):1241–1256
 21.
Boswell GP, Jacobs H, Ritz K, Gadd GM, Davidson FA (2007) The development of fungal networks in complex environments. B Math Biol 69(2):605–634
 22.
Coradin JH, Braun A, Viccini G, Krieger N, Mitchell DA (2011) A threedimensional discrete latticebased system for modeling the growth of aerial hyphae of filamentous fungi on solid surfaces: a tool for investigating microscale phenomena in solidstate fermentation. Biochem Eng J 54(3):164–171
 23.
Boswell GP, Davidson FA (2012) Modelling hyphal networks. Fungal Biol Rev 26(1):30–38
 24.
Heaton L, Obara B, Grau V, Jones N, Nakagaki T, Boddy L et al (2012) Analysis of fungal networks. Fungal Biol Rev 26(1):12–29
 25.
Nieminen L, Webb S, Smith MCM, Hoskisson PA (2013) A flexible mathematical model platform for studying branching networks: experimentally validated using the model actinomycete. Streptomyces coelicolor. PloS one 8(2):54316
 26.
Nielsen J (1992) Modelling the growth of filamentous fungi. In: Fiechter A (ed) Adv Biochem Eng Biotechnol, vol 46. Springer, pp 187–223
 27.
BartnickiGarcia S, Hergert F, Gierz G (1989) Computer simulation of fungal morphogenesis and the mathematical basis for hyphal (tip) growth. Protoplasma 153(1–2):46–57
 28.
Gierz G, BartnickiGarcia S (2001) A threedimensional model of fungal morphogenesis based on the vesicle supply center concept. J Theor Biol 208(2):151–164
 29.
Tindemans SH, Kern N, Mulder BM (2006) The diffusive vesicle supply center model for tip growth in fungal hyphae. J Theor Biol 238(4):937–948
 30.
Goriely A, Tabor M (2003) Biomechanical models of hyphal growth in actinomycetes. J Theor Biol 222(2):211–218
 31.
Boudaoud A (2003) Growth of walled cells: From shells to vesicles. Phys Rev Lett 91:018104
 32.
CaballeroLima D, Kaneva IN, Watton SP, Sudbery PE, Craven CJ (2013) The spatial distribution of the exocyst and actin cortical patches is sufficient to organize hyphal tip growth. Eukaryot cell 12(7):998–1008
 33.
Sugden KEP, Evans MR, Poon WCK, Read ND (2007) Model of hyphal tip growth involving microtubulebased transport. Phys Rev E 75(3):031909
 34.
Heaton LM, López E, Maini PK, Fricker MD, Jones NS (2010) Growthinduced mass flows in fungal networks. Proc R Soc Lond Ser B 277:3265–3274
 35.
Agrawal D, Bernstein P, Bertino E, Davidson S, Dayal U, Franklin M et al (2012) Challenges and opportunities with big data. CRA white paper
 36.
Harris SD, Read ND, Roberson RW, Shaw B, Seiler S, Plamann M et al (2005) Polarisome meets spitzenkörper: microscopy, genetics, and genomics converge. Eukaryot Cell 4(2):225–229
 37.
Riquelme M (2013) Tip growth in filamentous fungi: a road trip to the apex. Ann Rev Microbiol 67:587–609
 38.
Lew RR (2005) Mass flow and pressuredriven hyphal extension in Neurospora crassa. Microbiology 151(8):2685–2692
 39.
Lew RR (2011) How does a hypha grow? The biophysics of pressurized growth in fungi. Nat Rev Microbiol 9(7):509–518
 40.
Bizukojc M, Ledakowicz S (2006) A kinetic model to predict biomass content for Aspergillus niger germinating spores in the submerged culture. Process Biochem 41(5):1063–1071
 41.
Fiddy C, Trinci APJ (1976) Nuclei, septation, branching and growth of Geotrichum candidum. J Gen Microbiol 97:185–192
 42.
MouriñoPérez RR, Roberson RW, BartnickiGarcía S (2006) Microtubule dynamics and organization during hyphal growth and branching in Neurospora crassa. Fungal Genet Biol 43(6):389–400
 43.
RamosGarcía SL, Roberson RW, Freitag M, BartnickiGarcía S, MouriñoPérez RR (2009) Cytoplasmic bulk flow propels nuclei in mature hyphae of Neurospora crassa. Eukaryot Cell 8(12):1880–1890
 44.
Schunck T, Herrero S, Fischer R (2011) The Aspergillus nidulans cenpe kinesin kipa is able to dimerize and to move processively along microtubules. Curr Genet 57(5):335–341
Acknowledgements
The author would like to thank Frederik Spanhoff, Arthur Ram and Vera Meyer for supplying unpublished vesicle measurements and VM who inspired this work.
Compliance with ethical guidelines
Competing interests The author declares that he has no competing interests.
Author information
Appendices
Appendix A: Gradient of the radius
From Eq. 3,
Implicit differentiation leads to
Appendix B: Derivation of the generic model
Expanding the last term in Eq. 9 into a Taylor series results in
This yields upon substitution and neglect of terms in \((dx)^n\), \(n \ge 2\)
In the distal part, the balancing volume \(V(x_d)=V_d\) stays constant, i.e., \(\dot{V}(x_d) = \dot{V}_d = 0\). In the tip region, however, a stationary observer sees a growing volume \(\dot{V}(x_t) = \dot{V}_t\) as the apex moves away and the local radius increases. Hence, generally, with \(m_{\mathcal{S}_i} = V_i \mathcal{S}_i\)
or
The surface area \(A_i\) and volume \(V_i\) will differ depending on the part considered.
Taylor series of the volume balance, see Eq. 11, which is only affected by the cytoplasmic flow, leads to
Rearrangement gives
and
Combining it with Eq. 29 leads to the generic model
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Cite this article
King, R. A framework for an organellebased mathematical modeling of hyphae. Fungal Biol Biotechnol 2, 5 (2015) doi:10.1186/s4069401500142
Received:
Accepted:
Published:
Keywords
 Morphological model
 Hypha
 Vesicles