Vesicle transport and growth dynamics in Aspergillus niger : Microscale modeling of secretory vesicle flow and centerline extraction from confocal fluorescent data

In this paper, we present a mathematical model to describe filamentous fungal growth based on intracellular secretory vesicles (SVs), which transport cell wall components to the hyphal tip. Vesicular transport inside elongating hyphae is modeled as an advection – diffusion – reaction equation with a moving boundary, transformed into fixed coordinates, and discretized using a high ‐ order weighted essentially nonoscillatory discretization scheme. The model describes the production and the consumption of SVs with kinetic functions. Simulations are subsequently compared against distributions of SVs visualized by enhanced green fluorescent protein in young Aspergillus niger hyphae after germination. Intensity profile data are obtained using an algorithm scripted in ImageJ that extracts mean intensity distributions from 3D time ‐ lapse confocal measurement data. Simulated length growth is in good agreement with the experimental data. Our simulations further show that a decrease of effective vesicle transport velocity towards the tip can explain the observed tip accumulation of SVs.

extend at the tip and proteins are secreted into the surrounding medium (Berepiki, Lichius, & Read, 2011;Lanzetti, 2007;Riquelme, 2013;Steinberg, 2007;Virag & Harris, 2006). The SV fusion with the cell membrane is assumed to be mediated by the soluble NSF attachment protein receptor where NSF stands for N-ethyl-maleimide-sensitive fusion protein (SNAREs; Chen & Scheller, 2001). For this purpose, v-SNARE proteins are included in the vesicle membrane, while the corresponding t-SNAREs are placed into the cell membrane as the target membrane.
The SVs in A. niger can be tracked using the reporter strain FG7, which expresses an enhanced green fluorescent protein (eGFP) fused with the v-SNARE protein SncA located in the vesicle membrane enclosing the vesicle (eGFP-SncA; Kwon et al., 2014).
Spore germination is a crucial phase in any industrial fermentation. Tip extension highly relies on a coordinated transport to and from the hyphal apex, with SVs as the main tools of cargo delivery. In this study, the temporal change of SV concentration in the early growth phase after germination of A. niger is measured by image analysis using the algorithm presented. For interpreting the SV data, we developed a simplistic advection-diffusion-reaction model for SV transport to test whether deceleration of SVs towards the hyphal tip can be used to explain the observed SV accumulation at the tip.

| Modeling SV transport
Approaches to modeling cargo delivery to the hyphal tip vary immensely in the scale of the specific cell compartments covered.
They can be categorized according to the level of detail included, namely as microscale models for phenomena occurring in a single hypha or hyphal tip, mesoscale models describing hyphal network transport and macroscale models for entire fermentation processes (Sugai-Guérios, . Mechanistic models that couple transport and growth dynamics are typically located in between micro-and mesoscale model types. Only a few studies have specifically simulated production, transport and consumption of SVs inside a single, unbranched hypha. Prosser and Trinci (1979), as one of the earliest contributions, describe hyphal growth on an agar surface by vesicle production and convective transport to the tip. The spore is described as a simple cell compartment with a SV production term. Branching is considered as a result of vesicle accumulation when the flow is stopped behind the septum. The local resolution is quite coarse, and first-order finite differences are used. Furthermore, the predicted SVs distributions are not confirmed by measured values. Yang, King, Reichl, and Gilles (1992) calculate tip extension for a filamentous bacterium based on the balance of a limiting key component produced along the length of the hypha and consumed at the tip. The material flow of the key component is treated as being solely due to diffusion. Again, the focus is more on the development of mycelium network structures than on the single hypha. Although the model does well reflect the length growth of the primary branch, it lacks a real, biological counterpart for the key component, which, therefore, was not measured. López-Isunza, Larralde-Corona, and Viniegra-González (1997) developed a holistic model to simulate convective and diffusive transport of wall precursors and substrate in the cytosol of single A. niger hyphae growing on solid medium. They further included substrate diffusion through the solid medium and subsequent uptake following Michaelis-Menten kinetics. Precursor material is formed intracellularly from the substrate and depleted due to the production of cell wall material, that is for elongation. For the simulation of a young germling, the model would have to be extended to include the fungal spore at its boundary. The substrate and precursor concentrations simulated were not verified by measurements; thus, they must be regarded as hypothetical at that point. The SV concentration in the tip region of A. niger hyphae is typically increased. Assuming SVs to be the primary tool of cargo delivery, this should, as well, apply for precursor material, but was not replicated in their simulation.
A later model by Balmant et al. (2015) refines their model for single, unbranched aerial hyphae by balancing intracellular maltose concentrations in tank-like compartments along the length of the hypha. They discriminate between a distal nonvesicle-producing zone, a vesicle-producing zone and the transport towards the tip, where vesicles are consumed for tip extension. Production and consumption zones are kept at a constant length, following experimental observations by Trinci (1971). Neither SV nor nutrient concentrations were fit to measured data of A. niger. As the model describes the growth of developed aerial hyphae, the tip compartment alone measures at least 10 µm. Thus, the early developmental stages, where the first outgrowth appears and the hyphal tip is yet to be formed, can hardly be described.
Although the longitudinal growth behavior of filamentous fungi has been simulated in the publications mentioned above, none of the studies modeled the early developmental phase of the tip or the formation of the SV tip cluster, which is typical for actively growing hyphae. To fill this gap, our work proposes a basic, extendable mathematical model and compares it to length and SV measurements of an A. niger hypha.
This study is structured as follows: First, an image processing algorithm to extract local fluorescence distributions along the centerline of the fungal hypha is described. A quantitative mathematical description of the early germination period in filamentous fungi is then derived and the numerical framework to handle the resulting moving boundary system is presented. Subsequently, exemplary data of the intracellular movement of SVs is processed accordingly and compared against simulation results. Finally, we study the sensitivities of the parameters and initial conditions on length growth and on the SV concentration profiles of the simulation. and Ram (2010). Cultivation was done in μ-slides VI 0.4 , precoated with poly-L-lysine (ibidi GmbH, catalogue 80604) to ensure spore attachment on the bottom surface of the slide. The inoculated slides were incubated at 37°C, both during 3 hr swelling time and subsequent image acquisition.
The 3D time-series data stacks of germinating spores of A. niger were acquired via a confocal laser scanning microscope (CLSM; Leica SP8), utilizing an OPSL laser at λ = 488 nm for fluorophore excitation and a detector in the range of 490-547 nm to capture emission from eGFP-tagged vesicles (green light emission). Autofocus settings were developed to automatically capture the spores and to stack images equidistantly in a fixed range around them. A laser at λ = 552 nm wave length, 695-795 nm detector range, was used to avoid extensive photo bleaching during the autofocus procedure, taking advantage of the autofluorescence of the cell wall and avoiding fluorophore excitation. Images were taken with a water immersion objective HC PL APO CS2 63 × /1.20. The voxel depth was set to d Vxl = 1.8 µm, while the pixel width/length was set to w Pxl = 0.1 µm. A lower z-resolution was accepted as a compromise between spatial and temporal resolution, while maintaining an acceptable signal-to-noise ratio.

| Image analysis
The mathematical model ought to be compared against experimental data of the temporal change of relative fluorescent signal intensities of eGFP::SncA in A. niger. The data obtained from confocal laser scanning microscopy measurements in this study consists of threedimensional (3D) image stacks at distinct time points, with one or more fluorescent channels detected. For fitting the model, these data sets are required to be reduced into one dimension by averaging fluorescent image stacks along the discrete medial axis of the hyphae.
The challenge of extracting a centerline out of blurry, noisecorrupted confocal image data is interspersed with pitfalls of different kinds. It has been tackled in different segmentation algorithms in the past, most of which originate from medical imaging of human cells, particularly of blood vessels, neurovascular structures or neurites. Referring to the review of Kirbas and Quek (2004), who distinguish between six algorithmic classes, the extraction algorithm developed in this contribution belongs to the first category of classic pattern recognition techniques. It combines elemental filter and segmentation steps in Fiji (Schindelin et al., 2012) with a minimal amount of user interaction necessary.
The CLSM image stacks are preprocessed by background subtraction, denoising by median-filtering and iterative deconvolution.
Isotropy is set up by a quadrature interpolation-based reslicing procedure. Furthermore, drift in all spatial directions is compensated for by the realignment plugin "align slices in stack" (Rueden et al., 2017).
is again transformed according tõ̃̃̃̃̃̃̃x y z x y z , , , , Figure 1g. As a result, the hyphal tube is straightened in all spatial directions, with the image dimensions This final stack is used to extract the intensity profile along̃ỹ , averaging orthogonally iñx and̃z , Figure 1h. The steps for skeletonization are displayed in Figure 1i. The one-dimensional mathematical model presented later uses x as the spatial coordinate.
Therefore, the axial coordinate of a straightened hyphãỹ introduced above will be renamed as x in what follows to ease comparison.
The segmentation based on a thresholded projection is not unique, that is, cases can be constructed where either the point of greatest distance from the spore is not the tip but a tight inner curve perpendicular to the view or the longest skeleton branch cannot easily be identified from the projected image (e.g., a star-like structure). The algorithm might, in these cases, fit the centerline to a different centerline, shorter than the longest. Early tip growth after germination, as measured in our experiment, is, however, characterized by rather linear path trajectories, which reduces the risk of the above case occurring.

| Modeling SV tip accumulation
Transport inside a cellular environment can be driven by different mechanisms. Passive diffusion through the cytosol allows molecules to be transferred by thermally driven, random movements but essentially countervails an existing concentration gradient that can be observed in filamentous fungi. Active diffusion, as another mechanism, is caused, for example, by the activity of molecular motors that mix the cytoplasm (Brangwynne, Koenderink, MacKintosh, & Weitz, 2008). By contrast, directed active transport with molecular Centerline extraction algorithm based on confocal fluorescent image data (measured in CLSM), schematic view: (a) input: stacked time series of single conidium with germ tube, prefiltered, (b) sum projection (z-wise) of Nth time frame, (c) centerline extraction according to skeletonization pipeline (see (h)), (d) straightening along the skeleton in xy plane, (e) sum of projection in w, (f) centerline extraction, (g) straightening along the skeleton in the yz plane, (h) extraction of mean intensities oñ̃̃xz planes orthogonal tõỹ results in mean centerline intensity as a function of the hyphal length, and (i) skeletonization pipeline: thresholding, radius estimation with tubeness ImageJ plug-in (Rueden et al., 2017), skeletonization and centerline extraction [Color figure can be viewed at wileyonlinelibrary.com] motor proteins (such as the kinesins), along microtubules allows the cell to transfer cargo at the cost of adenosine triphosphate (ATP).
Furthermore, an unbound molecule can be transported passively inside the intracellular stream of the cytoplasm. In A. niger, active transport is the determining mechanism behind the long-distance transport and accumulation of SVs (Steinberg, 2007).
Individual SVs move at a pace above the velocity of hyphal extension. Near the tip, their concentration is not constant but shows a gradient along the length of the hypha, as apparent in microscopic images (see Figure 3a). We will, hereafter, assume that the transport speed on subapical microtubules exceeds the velocity expected on apical actin cables, that is the velocity reduces towards the tip. Crosslinking between both structures is very likely (Elie et al., 2015). Thus, the transition between both velocities will be rather smooth. In this article, we investigate whether a linear decrease in transport velocity towards the hyphal tip can be used to explain the formation of the observed SV tip accumulation.

| Model assumptions
Our dynamic model describes the active transport of SVs along microtubules (MT) inside the hyphal tube as well as diffusive vesicle motion at its tip. During the time period covered by the model, neither septation nor branch initiation has taken place and intracellular vesicle flow is neither stopped nor decreased elsewhere than in the tip region. Substrate uptake or secretion are not considered explicitly but can be included in the model structure. Furthermore, we do not distinguish between microvesicles and macrovesicles, which carry different kinds of enzymes, synthesizing the cell wall (Riquelme, 2013). Hence, all vesicles are assumed to be fully loaded (or at a constant ratio). Recycling of SVs after cargo release is not integrated into the model.
SVs will hereafter be represented as continuous, local concentrations x ( ) as the product of the local microtubule density F I G U R E 2 (a) Algorithmic details during the stretch phase. I: Young hyphal tip after germination, right tip boundary moves away; II: while the grid is stretched equidistantly, secretory vesicle fluxes, local production, and consumption are balanced inside a volume element dV; III: developed consumption zone Δx C and deceleration zone Δx u stay stationary relative to the hyphal tip. (b) Development of the tip velocity gradient during outgrowth: Velocity at the spore surface is held constant, while velocity at the tip u min (t) decreases to a minimum until condition L(t) = Δx u is met. Subsequently, the deceleration zone moves in parallel to the tip [Color figure can be viewed at wileyonlinelibrary.com] transport velocity at the very tip is assumed to be close to 0.
Therefore, a decrease in diameter would not show strong effects on the simulation. The radius is, thus, held constant mathematically at a mean measured radius of r = 3.73/2 μm.

| Balancing equations
Starting with a general mass balance for the infinitesimal small vo- results, with a vesicle mass flow ṁ entering the balance volume dV at x and leaving it at x dx + . Volumetric vesicle production and consumption terms are given by μ x t , P ( ) and μ x t , C ( ), respectively (see Figure 2a, II). Since SVs are produced inside the organism and cannot cross cell wall boundaries, there is no transport from or to the environment.
The mass flow of vesicles ṁ can be replaced by a flux density term in the x-direction, j x , or more refined by the vesicle concentration x t , ( ) by combining the advective transport with velocity u x t , ( ) due to active movement and the diffusion of unbound vesicles: x , .
Expansion of the second term on the right-hand side of Equation (6) into a Taylor series gives while neglecting all terms in (dx) n , n ≥ 2 leads to , .
Inserting m x t x t dV , , ( ) = ( )⋅ , differentiating with respect to t and dividing by dV A dx = ⋅ , yields The first term on the right-hand side describes the change of concentration due to internal fluxes, while the second captures source and sink terms due to vesicle production and consumption,

| Moving boundary problem
New balance volume is formed while the tip grows out of the spore ( Figure 2a), but existing SVs are already present inside the germ tube.
The problem is treated like a moving boundary problem to solve the transport equation on the growing balance domain (while keeping a confined number of numerical grid points and, therefore, equations to solve after subsequent discretization). Hence, the mathematical system is transformed according to Crank's front-fixing method (Crank, 1987), which artificially fixes the grid on the dimensionless axial coordinate η 0, 1 with the moving boundary at η 1 = .
The relationship η f x t f t , , ( ) = ( ) must be maintained for η x L = for an arbitrary function in the x-coordinate f x t , ( ) to be transformed into η f t , ( ). The basic model equation (9) is transformed accordingly into fixed coordinates (see the full derivation in Appendix B): Consequently, the total number of grid points N in the numerical grid stays constant when the hyphal tip extends. The transformative rule is, therefore, equivalent to an equidistant stretching of the numerical grid, while inserting an artificial flux term to cope with the increase in length. For simplicity, all following expressions concerning kinetics and boundary/initial conditions will be given as functions of the x-coordinate, but solved in this numerical scheme as functions of η.

| Initial and boundary conditions
The system is well-posed only if proper boundary and initial conditions are defined. As swelling conidia show considerable secretory vesicle activity, it can be reasonably assumed that vesicle flux from the spore contributes to tip growth in the very first period after germination. Since the vesicle flow from the spore into the germ tube is unknown and its contribution to tip growth decreases with increasing hyphal length, imposing a specific mass flux at the left boundary does not seem a sensible solution. Instead, we can set the boundary concentration v t 0, 1 ( ) = directly if we normalize the measured intensity profiles by the stationary left boundary value at later time steps, that is when the tip body has fully developed.
However, when the germ tube is still short, the left boundary value is enhanced by the evolving vesicle cluster, which is depleted when the tip moves forward after germination. As a first approach, we assume that the left boundary converges to v t 0, 1 ( ) = according to a simple modelv At the end of the germ tube, SVs are either fused with the cell membrane or recycled internally. Therefore, the right boundary is assumed to be perfectly absorbing, and the concentration is fixed at v 0 x L = = . The first measured hyphal length value and concentration profile are taken as initial length and initial distribution, respectively.

| Kinetics
For Equation (11), expressions for consumption, tip growth, production, transport and diffusion have to be specified.

Consumption
We assume that the hyphal extension rate links to the vesicle fusion at the tip membrane, that is faster growth necessitates higher tip consumption of SVs. The consumption term μ̃x t , C ′ ( ) describes this mechanism. Consumption by cell membrane fusion is further assumed to be restricted to a consumption zone x C Δ behind the apex where vesicles merge with the extending cell membrane (Figure 2a). It exhibits saturation kinetics because a carrier complex mediates SV fusion: where K C is the empirical saturation constant and k C is the maximal elongation velocity. We calculated Δ C = 13 µm by scaling the zone length obtained in Balmant et al. (2015) to the diameter of A. niger.

Tip growth
Elongation of the first hypha after germination typically undergoes a transition from exponential to linear length growth dynamics (Fiddy & Trinci, 1976;Gougouli & Koutsoumanis, 2013;King, 2015;Nanguy, Perrier-Cornet, Bensoussan, & Dantigny, 2010) and is modeled by a Monod-type kinetic expression, as a function of the sum of consumed SV in all grid cells within the consumption zone x C Δ . Furthermore, we add the flow of SVs that pass from the penultimate (i.e., the N 1 − th) node to the absorbing right boundary: where Y L is the extension of hyphal length per (dimensionless) mass of vesicle consumed.

Production
We expect vesicle synthesis to take place in the last x P Δ behind the tip of the young germ tube. Out of a general range of 5 m 2 0 m P µ ≤ Δ ≤ µ for filamentous fungi, as stated by Sugai-Guérios et al. (2015), we choose the maximum Δ P = 20 µm as an estimate. As a first attempt, we propose a fixed production rate, θ P , following Steinberg and Schuster (2011), who observed a relatively homogeneous distribution of Golgi equivalents (and, therefore, SV sources) in the tip area of A. nidulans:

Transport
As a first approach, transport speed is set to be at a constant value u max within the subapical hypha and to decrease linearly from this value to u t min ( ) at the very apex within the deceleration region x u Δ : During the early growth phase after germination, when x u Δ is not yet established, u t min ( ) decreases over time, until the condition L t x u ( ) = Δ is met: velocities would lead to a rapid accumulation in the first actin element in the simulation, which was deemed unlikely. Moreover, numerical issues due to the step in concentration are avoided by assuming a gradual reduction of the velocity.

Diffusion
If SVs were considered completely unbound and free to move inside the cytoplasm, the Einstein-Smoluchowski relationship would probably hold (Einstein, 1926(Einstein, /1956). However, studies revealed that vesicle-sized particles move significantly more slowly in the cytosolic environment (Luby-Phelps, 1999). Furthermore, as current scientific findings suggest an actin-mediated, star-shaped transport pattern in the apical dome, effective diffusion will merely serve as a simplified way to model isotropic cargo redistribution with the diffusion constant D = 5 µm 2 ·min −1 . Diffusion is expected in the whole compartment behind the spore, due to the lack of septation in the young germ tube (Trinci, 1974).

| Parameter identification procedure
The model parameters were identified by fitting the model to the measured SVs profiles and the length growth of the hypha. To ensure a priori identifiability of the model, some parameters were fixed either by approximate calibration or from literature data. This was necessary because several parameters (like k C and Y L ) were mathematically dependent. The parameters θ θ k x , , Δx C Width of consumption zone 13.00 µm Balmant et al. (2015) Y L Extension factor length growth 2 × 10 −2 µm·Ves −1 Fixed further. Parameter confidence intervals were computed using a Bayesian importance sampling approach (MCMC, viz., parallel tempering, using the PESTO toolkit by Stapor et al., 2018). As an objective function for parameter estimation, a combination of an unweighted least-squares formulation for vesicle concentration and length was calculated

| Model implementation details
The partial differential equation (11) is solved using the method of lines, that is by the discretization and calculation of spatial derivatives with finite difference stencils. As tip-accumulated vesicles typically show sharp gradients, standard linear finite difference stencils (such as upwind or central) will result in numerical artifacts and cannot be used. Therefore, high-order spatial derivatives are calculated using polynomial reconstruction with pre-computed weights, rooting to the work of Shu (1998). The resulting system of ODEs is solved using an explicit Runge-Kutta routine in C, which was compiled and executed in MATLAB. We used an Intel Core i7-6700 3.40 GHz CPU, 32 GB RAM system for all calculations involved in this study.

| RESULTS
Model parameters were fit to a germ tube that was representative of the population. The predicted length growth (simulated mean) of A.
niger is compared to measured values in Figure 3b. Uncertainty bounds for the simulation are plotted in shades of blue, based on the empirical quantiles of the predicted length growth and vesicle distributions outlined above. Simulative and measured length data points match closely. Length is overestimated slightly at times t < 120 min; however, all measured values lie within predicted uncertainty bounds at all times. This overestimation carries over to the comparison of measured and simulated SV evolution that shows an excellent fit of the space-dependent concentration profile, but lacks the correct length. Figure 3c shows The parameter estimates θ obtained from the identification procedure are shown in Table 1. The estimated length of the velocity deceleration zone x u Δ spans the last 4.85 µm of the germ tube. The estimated production factor is θ P = 1.10 Ves·µm −3 min −1 . The estimate for the maximum specific growth rate shows high uncertainty with a relative error of 61.0%. The MCMC results (Appendix C) prove a negative correlation between the production factor θ P and k Ĉ : When more SVs are produced, fewer SVs are consumed in order not to exceed elongation velocity. The initial SV profile was multiplied by a factor f IC to test the influence of varying initial conditions, as plotted in Figure 4a. The length deviation scales linearly, since more SVs linearly translate into KUNZ ET AL.

| Sensitivity analysis
| 2883 length growth, following Equation (14). The effect of f IC on hyphal length, however, is low compared to the effects of the model parameters.

| DISCUSSION
In this article, we derived a model to describe the build-up of the SV tip cluster in young A. niger hyphae after germination. The model was fit to intracellular SV concentration measurements, which were obtained from 3D confocal laser scanning microscopy data by using the centerline algorithm described. This processing algorithm developed extracts mean axial distributions from 3D time-lapse data of hyphal objects and was used to quantify SVs inside young A. niger hyphae.
The algorithm could be used to quantify the density of other cell organelles from fluorescent confocal image stacks and will set the data basis for future transport modeling approaches. However, the algorithm places relatively high demands on the position of the sample objects (in this case: of the individual spores and the hyphae), as they need to be clearly separated from each other in the image. In this experiment, the germinating spores formed relatively few clusters, which was due to the low spore density.
The mathematical model describes fungal tip extension and transport of SVs in young A. niger hyphae, which develop from spores or as new branches in growing mycelia. The model uses a fixed rate to describe SV production and saturation kinetics to model SV consumption due to fusion with the cell membrane. Both SV production and consumption terms are restricted to distinct zones behind the apex, which ensures that individual hyphal tips can grow independently and on self-supply. The growth equation in the model is a function of the mass of vesicles fusing with the cell membrane. As a first step, we wanted to test the hypothesis that a simple decrease in transport velocity towards the apex could be responsible for the increased SV occurrence in the apical dome. With the overall model, the key phenomenon of SV accumulation at the hyphal apex, as well as the speed of tip extension, could be replicated. The model showed weaknesses in the description of the early formation of the tip profile, see Figure 3c.
A mathematically similar approach has been used to model the growth of fungal aerial hyphae . Although the model has not been validated for A. niger and its tank-like discretization is too coarse for hyphae in such early stages of development, it would be interesting to compare its abilities to describe SV accumulation in fully grown hyphae.
The parameter estimation result proved the overall ability of the model to describe the dynamics of the establishment of high vesicle density at the tips of young germ tubes in A. niger. Our assumption that the effective transport velocity decreases towards the tip led to transport dynamics in good agreement with the SV measurement data. An extended number of young germ tubes would have to be taken into account to quantify the apparent growth and transport heterogeneity among outgrowing hyphae in future studies. Notably, uncalibrated confocal microscopy can only provide relative values of the vesicle concentration, which depend on multiple experimental factors. On the one hand, this allowed us to normalize the image data.
On the other, the kinetic parameters identified under these conditions are also normalized values and, therefore, lack the ability to predict absolute vesicle densities. In situ measurements of the absolute SV concentrations in the cytosol could deal with these restrictions. It will serve to quantify variation in SV distribution between static and flow cultivation conditions to address the ques- strategy when choosing the solver. All these issues will be dealt with in future studies.