TY - JOUR
AU - Boswell, Graeme, P
AB - Abstract Nearly all life forms require iron to survive and function. Microorganisms utilize a number of mechanisms to acquire iron including the production of siderophores, which are organic compounds that combine with ferric iron into forms that are easily absorbed by the microorganism. There has been significant experimental investigation into the role, distribution and function of siderophores in fungi but until now no predictive tools have been developed to qualify or quantify fungi-initiated siderophore–iron interactions. In this investigation, we construct the first mathematical models of siderophore function related to fungi. Initially, a set of partial differential equations are calibrated and integrated numerically to generate quantitative predictions on the spatio-temporal distributions of siderophores and related populations. This model is then reduced to a simpler set of equations that are solved algebraically giving rise to solutions that predict the distributions of siderophores and resultant compounds. These algebraic results require the calculation of zeros of cross products of Bessel functions and thus new algebraic expansions are derived for a variety of different cases that are in agreement with numerically computed values. The results of the modelling are consistent with experimental data while the analysis provides new quantitative predictions on the time scales involved between siderophore production and iron uptake along with how the total amount of iron acquired by the fungus depends on its environment. The implications to bio-technological applications are briefly discussed. mathematical model, partial differential equations, numerical solution, ferric iron uptake 1. Introduction Iron is an essential element for nearly all life forms. In humans, iron deficiency can lead to several chronic medical conditions (such as anemia, Zimmermann & Hurrell, 2007; Beard, 2008), whereas in plants, insufficient amounts of iron can severely hinder growth, which is particularly problematic since one third of the world’s soils are considered to be iron deficient due to the insolubility of ferric iron present in the environment (Marschner, 1995). Indeed, nutritional iron is not readily available in the terrestrial environment, and thus, microorganisms have evolved mechanisms to cope with its scarcity by developing processes to convert and subsequently uptake iron to aid in their growth. These mechanisms have been studied at the molecular level for various microscopic eukaryotes including bacteria and pathogenic fungi (Philpott et al., 2012). In fungi, four different mechanisms for the acquisition of iron have been identified (e.g. Van der Helm et al., 1994; Renshaw et al., 2002; Haas, 2014 and references therein): (i) shuttle mechanism: ferric iron uptake mediated by ferric iron-specific chelators (siderophores), (ii) direct-transfer mechanism: reductive iron assimilation, (iii) esterase-reductase mechanism: low-affinity ferrous iron uptake and (iv) reductive mechanism: heme uptake and degradation. In this work, we focus attention on the first, and most common, of these mechanisms. Under iron-limited conditions, many microorganisms produce and secrete small organic molecules called siderophores (Schwyn & Neilands, 1987; Saha et al., 2016). Siderophores are low molecular weight iron chelating compounds that move by Brownian motion and have a high affinity for ferric iron. Once the siderophores are attached to the ferric iron, the siderophore–iron complexes are transported by diffusion (Srivastava et al., 2013) and can be acquired by the organism, whereupon the iron is internalized and used to support further biomass growth and function. Siderophores have drawn much attention in recent times due to their potential roles and applications in various bio-technologies including agriculture, ecology, bio-remediation, bio-control, bio-sensor and medicine (Saha et al., 2016). Their significance in applications are mainly due to siderophores having the ability to bind to a variety of metals in addition to iron (Braud et al., 2009; Bellenger et al., 2013; Sasirekha & Srividya, 2016). For example, siderophores play a crucial role in mobilizing metals from metal contaminated soils (Ahmed & Holmström, 2014 and references therein). Additionally in bio-control, microorganisms that produce certain siderophores can take up iron from around their immediate vicinity and invade a competitor’s space in search for iron, which leads to the suppression of growth of several fungal pathogens (McLoughin et al., 1992; Verma et al., 2011). Siderophores are classified by the ligands (an ion, molecule or molecular group that binds to another chemical entity to form a larger complex) used to chelate the ferric iron that can be categorized as catecholates, hydroxamates and carboxylates (Winkelmann, 1991, 2002; Ahmed & Holmström, 2014). Fungi mostly produce siderophores that fall in the ‘hydroxamates’ category and most species of fungi make more than one type of siderophore, possibly to adapt to different environmental conditions (Renshaw et al., 2002; Perez-Meranda et al., 2007; Johnson, 2008). Thus, various assays have been developed to detect the different phenotypes of siderophores. While these assays are useful for identifying various siderophores, numerous assays would have to be formed independently to detect all possible forms of siderophores, of which there are more than 500 known distinct types (Boukhalfa et al., 2003; Kraemer et al., 2005). Schwyn & Neilands (1987) developed a universal siderophore detection assay using chrome azurol S (CAS) and hexadecyltrimethylammonium bromide (HDTMA) as visual indicators of the presence and function of siderophores. The CAS/HDTMA complexes tightly bond with ferric iron and become blue in colour. When a strong iron chelator, such as a siderophore, removes iron from the dye complex, the colour typically changes from blue to either orange, magenta or purple, depending on the exact assay (Bertrand et al., 2010). The toxicity induced by the HDTMA indicator can, in certain species, inhibit and even prevent the normal growth and function of the fungus (Schwyn & Neilands, 1987). Consequently, numerous later studies (e.g. Milagres et al., 1999) have been based around a split Petri dish where the HDTMA indicator is added to one semi-circular region but absent from the other half; such configurations have been successfully modelled by one of the authors (Choudhury, 2019). Despite their widespread existence, there has been relatively little attempt at the mathematical modelling of siderophores and their interaction with iron. In fungi, their mathematical treatment has typically been focussed on quantifying siderophore extent using simple ad-hoc approaches, such as measuring the physical distance of the colour change on a Petri dish or placing square paper underneath the Petri dish and recording the change in area over a time period (Machuca & Milagres, 2003; Bogumił et al., 2013; Ghosh et al., 2015; Andrews et al., 2016a,b). However, siderophores produced by bacteria have received more advanced mathematical treatment, typically using sets of differential equations (e.g. Eberi & Collinson, 2009; Niehus et al., 2017). Leventhal et al. (2019) developed the most insightful mathematical model by considering siderophores produced from a single non-moving and isolated bacteria cell and their subsequent interaction with iron in a marine environment to form siderophore–iron complexes and represented this process using a simple reaction–diffusion equation. Consequently, and given the sheer volume of applications involving fungi described above, it is timely that such a mathematical modelling exercise is performed that focuses on siderophore production involving an expanding fungal colony and thus significantly extending previous treatments of siderophore function. In this article, a set of partial differential equations (PDEs) is developed that model the growth of a fungal biomass in response to nutrients and which produces siderophores to acquire iron from the external environment. The models are less concerned with how the biomass subsequently uses the iron, rather the models predict the quantity of iron acquired by the biomass and how iron is distributed in the external environment as a result of siderophore interactions, and thus provides quantitative predictions related to the experimental protocols described above. A mathematical model is developed in Section 2 that simulates the growth of a mycelium, the production of siderophores and their resultant interaction with iron in a planar domain, representing typical experimental protocol corresponding to the growth of a fungus in a Petri dish. The effect of different concentrations of iron and external nutrients are investigated by solving the equations numerically. These simulations motivate the construction of a simplified set of equations, considered in Section 3, which focus on the siderophore dynamics. Algebraic solutions are constructed that describe the temporal evolution of the siderophore dynamics towards a steady-state distribution and are consistent with the numerical approach. These algebraic solutions make use of various asymptotic expansions applied to cross-products of Bessel functions and hence new results and methods are developed accordingly. The implications of the results and future work are discussed in Section 4. 2. Siderophore–iron interactions from an expanding biomass 2.1 Model equations Due to the dense network structure of a fungal mycelium, a continuum approach is used to model its growth in a planar setting, representing mycelial expansion in a Petri dish. The growth and function of a fungus in such settings has been previously modelled by Boswell et al. (2003) and expanded upon in a series of papers (e.g. Boswell et al., 2007; Choudhury et al., 2018 and references therein). In short, a fungal mycelium comprises a network of tubes, termed hyphae, that can branch, extend at their unbounded ends, fuse with other hyphae (anastomosis), acquire new growth material from the external environment (uptake) and redistribute that material through the network (translocation). For the purposes of modelling, the mycelium is assumed to comprise three variables representing active hyphae (denoted by |$\rho$| and corresponding to those hyphae involved in nutrient uptake, branching, anastomosis and translocation), inactive hyphae (denoted by |$\rho ^\prime$| corresponding to hyphae no longer involved in colony function but still remaining part of the mycelium) and hyphal tips (denoted by |$n$|⁠) representing the expanding ends of active hyphae. Briefly, hyphal tips move predominantly in a straight line but with some random variations that is modelled by an advective process directed away from hyphae coupled with a diffusive process representing the random reorientation. (This growth characteristic is a consequence of the delivery of vesicles from the Spitzenkörper to the hemi-ellipsoid-shaped apical tip, Riquelme et al., 2018). New hypha are therefore formed from the trail left behind a tip as it moves and thus the tip flux corresponds to the creation of hyphal biomass. Thus, the magnitude of the flux is a convenient approximation of the amount of new material created through the movement of hyphal tips. Tips are created through branching along existing active hyphae and are lost through anastomosis also with active hyphae. It is assumed that a single generic substrate is responsible for growth. This substrate exists in two forms: external to the mycelium (with density |$s_e$|⁠) and held within the mycelium (with density |$s_i$|⁠). The external substrate may represent combinations of carbon, nitrogen and trace metals other than iron while the internalized substrate additionally includes iron. Internally located substrate is translocated through the biomass structure by a combination of diffusion and active transport directed towards the hyphal tips, the latter of which has a metabolic cost and there is a further cost associated with the movement of hyphal tips. Consistent with experimental evidence, tip flux and branching rates increase with the internal substrate (Gruhn et al., 1992), and this resource is also used to uptake external substrate. It is assumed that the biomass is in an iron-depleted state, and thus, siderophores are being released throughout its extent. Consistent with the nutrient uptake process, it is assumed siderophore production can only arise in the presence of sufficient energy reserves, and in the absence of experimental evidence to the contrary, it is assumed that siderophore production is proportional to the internal substrate concentration and the density of active biomass with |$r_1$| denoting the constant of proportionality. When released to the external environment, siderophores (denoted by |$C$|⁠) diffuse and bind with iron (denoted by |$I$|⁠) to form siderophore–iron complexes (denoted by |$V$|⁠) and standard enzyme-reaction kinetics are assumed to describe this binding process with |$r_2$| denoting the rate constant. These complexes subsequently diffuse and are absorbed by the biomass across hyphal cell walls. As previously mentioned, there are in excess of 500 different types of siderophores with quantitatively and qualitatively different characteristics, and consequently, there are a multitude of different pathways via which the fungus acquires iron from the siderophore–iron complexes (Howard, 1999; Winkelmann, 2007). Simple diffusion across the hyphal cell wall is clearly common to all, and hence, this process is used to account for the iron uptake, where |$r_3$| is the uptake rate constant. Once internalized, the iron forms a component of the generic internal substrate that is subsequently used to promote further growth via hyphal tip extension and translocation and to acquire additional resources, including more iron. The uptake and subsequent conversion of the siderophore–iron complex into the generic internalized substrate has an associated cost and hence the effective acquisition rate of the complex, |$r_3^\prime$|⁠, is less than the overall uptake rate, |$r_3$|⁠. Thus, the entire system can be modelled using the mixed hyperbolic-parabolic set of partial differential equations given by $$\begin{equation}\hskip-135pt \rho_t = \left| v s_i n \nabla \rho + D_n s_i \nabla n \right| - d_\rho \rho, \end{equation}$$(2.1a) $$\begin{equation}\hskip-198pt \rho_t^\prime = d_\rho \rho - d_i \rho^\prime, \end{equation}$$(2.1b) $$\begin{equation}\hskip-85pt n_t = \nabla \cdot \left( v s_i n \nabla \rho + D_n s_i \nabla n \right) + \alpha s_i \rho - \beta n \rho, \end{equation}$$(2.1c) $$\begin{align} s_{i_t} =&\ \nabla \cdot \left( D_i \rho \nabla s_i - D_a \rho s_i \nabla n\right) + c_1 \rho s_i s_e - c_2 \left| v s_i n \nabla \rho + D_n s_i \nabla n \right| \nonumber \\ &- c_4 \left| D_a \rho s_i \nabla n \right| - r_1 \rho s_i + r_3^\prime V \rho, \end{align}$$(2.1d) $$\begin{equation}\hskip-172pt s_{e_t} = D_e \nabla^2 s_e - c_3 \rho s_i s_e, \end{equation}$$(2.1e) $$\begin{equation}\hskip-182pt I_t = D_I \nabla^2 I - r_2 I C, \end{equation}$$(2.1f) $$\begin{equation}\hskip-148pt C_t = D_C \nabla^2 C + r_1^\prime \rho s_i - r_2 I C, \end{equation}$$(2.1g) $$\begin{equation}\hskip-148pt V_t = D_{V} \nabla^2 V + r_2 IC - r_3 V \rho. \end{equation}$$(2.1h) The model variables and parameters along with their calibrated values are given in Tables 1 and 2, respectively. Table 1 Summary of model variables used in equation (2.1). Variable . Description . Unit . |$\rho$| Active hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$\rho ^\prime$| Inactive hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$n$| Hyphal tip density tips cm|$^{-2}$| |$s_i$| Internal substrate concentration mol cm|$^{-2}$| |$s_e$| External substrate concentration mol cm|$^{-2}$| |$I$| Concentration of free iron mol cm|$^{-2}$| |$C$| Concentration of siderophores mol cm|$^{-2}$| |$V$| Concentration of siderophore–iron complex mol cm|$^{-2}$| Variable . Description . Unit . |$\rho$| Active hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$\rho ^\prime$| Inactive hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$n$| Hyphal tip density tips cm|$^{-2}$| |$s_i$| Internal substrate concentration mol cm|$^{-2}$| |$s_e$| External substrate concentration mol cm|$^{-2}$| |$I$| Concentration of free iron mol cm|$^{-2}$| |$C$| Concentration of siderophores mol cm|$^{-2}$| |$V$| Concentration of siderophore–iron complex mol cm|$^{-2}$| Open in new tab Table 1 Summary of model variables used in equation (2.1). Variable . Description . Unit . |$\rho$| Active hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$\rho ^\prime$| Inactive hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$n$| Hyphal tip density tips cm|$^{-2}$| |$s_i$| Internal substrate concentration mol cm|$^{-2}$| |$s_e$| External substrate concentration mol cm|$^{-2}$| |$I$| Concentration of free iron mol cm|$^{-2}$| |$C$| Concentration of siderophores mol cm|$^{-2}$| |$V$| Concentration of siderophore–iron complex mol cm|$^{-2}$| Variable . Description . Unit . |$\rho$| Active hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$\rho ^\prime$| Inactive hyphal density cm|$^{-1}$| (cm hyphae cm|$^{-2}$|⁠) |$n$| Hyphal tip density tips cm|$^{-2}$| |$s_i$| Internal substrate concentration mol cm|$^{-2}$| |$s_e$| External substrate concentration mol cm|$^{-2}$| |$I$| Concentration of free iron mol cm|$^{-2}$| |$C$| Concentration of siderophores mol cm|$^{-2}$| |$V$| Concentration of siderophore–iron complex mol cm|$^{-2}$| Open in new tab Table 2 Parameter values used in model equations (2.1) with initial data (2.2). The values are taken from 1Boswell et al. (2002); 2Boswell et al. (2003); 3Perez-Meranda et al. (2007); 4Eberi & Collinson (2009); 5Leventhal et al. (2019), while the remaining parameters were assumed to take values consistent with those in similar processes. The value of |$R_{dish}$| was chosen to represent a Petri dish of radius 2 cm for computational convenience. Parameter . Value . Description . Unit . |$v$| 0.5 Tip velocity|$^{2}$| cm|$^5$| day|$^{-1}$| mol|$^{-1}$| |$D_n$| 0.1 Tip diffusion|$^2$| cm|$^4$| day|$^{-1}$| mol|$^{-1}$| |$d_\rho$| 0.2 Hypha inactivation rate|$^1$| day|$^{-1}$| |$d_i$| 0 Inactive hypha decay rate|$^1$| day|$^{-1}$| |$\alpha$| 10 000 Branching rate|$^2$| cm mol|$^{-1}$| day|$^{-1}$| |$\beta$| 10 000 Anastomosis rate|$^2$| cm day|$^{-1}$| |$D_i$| 10 Internal substrate diffusion coefficient|$^2$| cm|$^3$| day|$^{-1}$| |$D_a$| 10 Internal substrate active transport|$^2$| cm|$^5$| day|$^{-1}$| |$c_1$| 900 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_2$| 1 Tip extension costs|$^1$| mol cm|$^{-1}$| |$c_3$| 1000 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_4$| |$10^{-8}$| Active translocation costs|$^{2}$| cm|$^{-1}$| |$D_e$| 0.0001 External substrate diffusion coefficient|$^{1}$| cm|$^2$| day|$^{-1}$| |$D_I$| 0.000864 Iron diffusion coefficient|$^4$| cm|$^2$| day|$^{-1}$| |$D_C$| 0.3 Siderophore diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$D_V$| 0.3 Complex diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$r_1$| |$10^{-7}$| Siderophore production costs|$^5$| cm day|$^{-1}$| |$r_1^\prime$| 100 Production of siderophores|$^5$| cm day|$^{-1}$| |$r_2$| 100 Complex production rate|$^5$| cm|$^2$| mol|$^{-1}$| day|$^{-1}$| |$r_3$| 1000 Complex uptake rate cm day|$^{-1}$| |$r_3^\prime$| 900 Conversion of iron to substrate cm day|$^{-1}$| |$R_{dish}$| 2 Radius of Petri dish cm |$R_{plug}$| 0.2 Radius of inoculum|$^1$| cm |$\rho _0$| 0.1 Initial biomass density|$^1$| cm|$^{-1}$| |$n_0$| 0.1 Initial tip density|$^1$| cm|$^{-2}$| |$s_{i_0}$| 0.4 Initial internal substrate density |$^{1}$| mol cm|$^{-2}$| |$s_{e_0}$| 0.6 Initial external substrate density |$^1$| mol cm|$^{-2}$| |$I_0$| 0.004 Initial iron concentration|$^{3}$| mol cm|$^{-2}$| Parameter . Value . Description . Unit . |$v$| 0.5 Tip velocity|$^{2}$| cm|$^5$| day|$^{-1}$| mol|$^{-1}$| |$D_n$| 0.1 Tip diffusion|$^2$| cm|$^4$| day|$^{-1}$| mol|$^{-1}$| |$d_\rho$| 0.2 Hypha inactivation rate|$^1$| day|$^{-1}$| |$d_i$| 0 Inactive hypha decay rate|$^1$| day|$^{-1}$| |$\alpha$| 10 000 Branching rate|$^2$| cm mol|$^{-1}$| day|$^{-1}$| |$\beta$| 10 000 Anastomosis rate|$^2$| cm day|$^{-1}$| |$D_i$| 10 Internal substrate diffusion coefficient|$^2$| cm|$^3$| day|$^{-1}$| |$D_a$| 10 Internal substrate active transport|$^2$| cm|$^5$| day|$^{-1}$| |$c_1$| 900 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_2$| 1 Tip extension costs|$^1$| mol cm|$^{-1}$| |$c_3$| 1000 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_4$| |$10^{-8}$| Active translocation costs|$^{2}$| cm|$^{-1}$| |$D_e$| 0.0001 External substrate diffusion coefficient|$^{1}$| cm|$^2$| day|$^{-1}$| |$D_I$| 0.000864 Iron diffusion coefficient|$^4$| cm|$^2$| day|$^{-1}$| |$D_C$| 0.3 Siderophore diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$D_V$| 0.3 Complex diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$r_1$| |$10^{-7}$| Siderophore production costs|$^5$| cm day|$^{-1}$| |$r_1^\prime$| 100 Production of siderophores|$^5$| cm day|$^{-1}$| |$r_2$| 100 Complex production rate|$^5$| cm|$^2$| mol|$^{-1}$| day|$^{-1}$| |$r_3$| 1000 Complex uptake rate cm day|$^{-1}$| |$r_3^\prime$| 900 Conversion of iron to substrate cm day|$^{-1}$| |$R_{dish}$| 2 Radius of Petri dish cm |$R_{plug}$| 0.2 Radius of inoculum|$^1$| cm |$\rho _0$| 0.1 Initial biomass density|$^1$| cm|$^{-1}$| |$n_0$| 0.1 Initial tip density|$^1$| cm|$^{-2}$| |$s_{i_0}$| 0.4 Initial internal substrate density |$^{1}$| mol cm|$^{-2}$| |$s_{e_0}$| 0.6 Initial external substrate density |$^1$| mol cm|$^{-2}$| |$I_0$| 0.004 Initial iron concentration|$^{3}$| mol cm|$^{-2}$| Open in new tab Table 2 Parameter values used in model equations (2.1) with initial data (2.2). The values are taken from 1Boswell et al. (2002); 2Boswell et al. (2003); 3Perez-Meranda et al. (2007); 4Eberi & Collinson (2009); 5Leventhal et al. (2019), while the remaining parameters were assumed to take values consistent with those in similar processes. The value of |$R_{dish}$| was chosen to represent a Petri dish of radius 2 cm for computational convenience. Parameter . Value . Description . Unit . |$v$| 0.5 Tip velocity|$^{2}$| cm|$^5$| day|$^{-1}$| mol|$^{-1}$| |$D_n$| 0.1 Tip diffusion|$^2$| cm|$^4$| day|$^{-1}$| mol|$^{-1}$| |$d_\rho$| 0.2 Hypha inactivation rate|$^1$| day|$^{-1}$| |$d_i$| 0 Inactive hypha decay rate|$^1$| day|$^{-1}$| |$\alpha$| 10 000 Branching rate|$^2$| cm mol|$^{-1}$| day|$^{-1}$| |$\beta$| 10 000 Anastomosis rate|$^2$| cm day|$^{-1}$| |$D_i$| 10 Internal substrate diffusion coefficient|$^2$| cm|$^3$| day|$^{-1}$| |$D_a$| 10 Internal substrate active transport|$^2$| cm|$^5$| day|$^{-1}$| |$c_1$| 900 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_2$| 1 Tip extension costs|$^1$| mol cm|$^{-1}$| |$c_3$| 1000 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_4$| |$10^{-8}$| Active translocation costs|$^{2}$| cm|$^{-1}$| |$D_e$| 0.0001 External substrate diffusion coefficient|$^{1}$| cm|$^2$| day|$^{-1}$| |$D_I$| 0.000864 Iron diffusion coefficient|$^4$| cm|$^2$| day|$^{-1}$| |$D_C$| 0.3 Siderophore diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$D_V$| 0.3 Complex diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$r_1$| |$10^{-7}$| Siderophore production costs|$^5$| cm day|$^{-1}$| |$r_1^\prime$| 100 Production of siderophores|$^5$| cm day|$^{-1}$| |$r_2$| 100 Complex production rate|$^5$| cm|$^2$| mol|$^{-1}$| day|$^{-1}$| |$r_3$| 1000 Complex uptake rate cm day|$^{-1}$| |$r_3^\prime$| 900 Conversion of iron to substrate cm day|$^{-1}$| |$R_{dish}$| 2 Radius of Petri dish cm |$R_{plug}$| 0.2 Radius of inoculum|$^1$| cm |$\rho _0$| 0.1 Initial biomass density|$^1$| cm|$^{-1}$| |$n_0$| 0.1 Initial tip density|$^1$| cm|$^{-2}$| |$s_{i_0}$| 0.4 Initial internal substrate density |$^{1}$| mol cm|$^{-2}$| |$s_{e_0}$| 0.6 Initial external substrate density |$^1$| mol cm|$^{-2}$| |$I_0$| 0.004 Initial iron concentration|$^{3}$| mol cm|$^{-2}$| Parameter . Value . Description . Unit . |$v$| 0.5 Tip velocity|$^{2}$| cm|$^5$| day|$^{-1}$| mol|$^{-1}$| |$D_n$| 0.1 Tip diffusion|$^2$| cm|$^4$| day|$^{-1}$| mol|$^{-1}$| |$d_\rho$| 0.2 Hypha inactivation rate|$^1$| day|$^{-1}$| |$d_i$| 0 Inactive hypha decay rate|$^1$| day|$^{-1}$| |$\alpha$| 10 000 Branching rate|$^2$| cm mol|$^{-1}$| day|$^{-1}$| |$\beta$| 10 000 Anastomosis rate|$^2$| cm day|$^{-1}$| |$D_i$| 10 Internal substrate diffusion coefficient|$^2$| cm|$^3$| day|$^{-1}$| |$D_a$| 10 Internal substrate active transport|$^2$| cm|$^5$| day|$^{-1}$| |$c_1$| 900 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_2$| 1 Tip extension costs|$^1$| mol cm|$^{-1}$| |$c_3$| 1000 Nutrient uptake rate|$^1$| cm|$^3$| mol|$^{-1}$| day|$^{-1}$| |$c_4$| |$10^{-8}$| Active translocation costs|$^{2}$| cm|$^{-1}$| |$D_e$| 0.0001 External substrate diffusion coefficient|$^{1}$| cm|$^2$| day|$^{-1}$| |$D_I$| 0.000864 Iron diffusion coefficient|$^4$| cm|$^2$| day|$^{-1}$| |$D_C$| 0.3 Siderophore diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$D_V$| 0.3 Complex diffusion coefficient|$^5$| cm|$^2$| day|$^{-1}$| |$r_1$| |$10^{-7}$| Siderophore production costs|$^5$| cm day|$^{-1}$| |$r_1^\prime$| 100 Production of siderophores|$^5$| cm day|$^{-1}$| |$r_2$| 100 Complex production rate|$^5$| cm|$^2$| mol|$^{-1}$| day|$^{-1}$| |$r_3$| 1000 Complex uptake rate cm day|$^{-1}$| |$r_3^\prime$| 900 Conversion of iron to substrate cm day|$^{-1}$| |$R_{dish}$| 2 Radius of Petri dish cm |$R_{plug}$| 0.2 Radius of inoculum|$^1$| cm |$\rho _0$| 0.1 Initial biomass density|$^1$| cm|$^{-1}$| |$n_0$| 0.1 Initial tip density|$^1$| cm|$^{-2}$| |$s_{i_0}$| 0.4 Initial internal substrate density |$^{1}$| mol cm|$^{-2}$| |$s_{e_0}$| 0.6 Initial external substrate density |$^1$| mol cm|$^{-2}$| |$I_0$| 0.004 Initial iron concentration|$^{3}$| mol cm|$^{-2}$| Open in new tab The flux term in equation (2.1c) corresponds to the motion of hyphal tips accounting for their straight line growth habit (where |$s_i$| accounts for the role of the growth promoting substrate in the process) coupled with variations about this orientation, modelled using diffusion. The parameter |$v$|⁠, corresponding to the straight line growth habit of individual hyphae, is influenced by toxicity in the external environment; in particular, tip extension can be inhibited through the presence of the HDTMA visual indicator used to detect the presence of siderophores (Schwyn & Neilands, 1987). Indeed, numerous studies (e.g. Fomina et al., 2000) have shown that the ability of fungi to colonize space occupied by toxic material is increased through the availability of nutrients such as carbon. Consequently, it is tacitly assumed that the HDTMA indicator is uniformly distributed and at a concentration that does not prevent the biomass from expanding so that |$v$| may be regarded as a positive constant and thus the expansion of the model biomass into the space where the HDTMA visual indicator is present is consistent with experimental observations. Furthermore, this phenomenon further justifies the explicit modelling of both an external substrate, representing nutrients that assist the fungi in overcoming the toxicity and the iron distribution. The metabolic cost of tip movement is accounted for in equation (2.1d) through the parameter |$c_2$|⁠, while the trail left behind the tip, and thus the creation of new hyphae, is given by the related term in equation (2.1a). The flux in equation (2.1d) represents movement of internally held material through the network (i.e. translocation) having both diffusive and directed components, the latter towards hyphal tips and having a metabolic cost. Equations (2.1a)–(2.1e) are precisely those in Boswell et al. (2003). In equations (2.1f)–(2.1h) the iron, siderophore and the siderophore–iron complex populations are assumed to undergo standard Fickian diffusion with coefficients |$D_I, D_C$| and |$D_V$|⁠, respectively. Note that the key function of siderophores is to increase the mobility of iron, which is achieved through the formation of siderophore–iron complexes. Thus |$D_I<D_V$|⁠. See Howard (1999) and Leventhal et al. (2019) for further details and discussion of calibration. Equations (2.1) are considered on a domain representing a typical experimental protocol, i.e. a circular Petri dish of radius |$R_{dish}$| with an initially uniform growth medium inoculated at its centre by a small circular plug of biomass of radius |$R_{plug}$|⁠. Consequently, the biomass is initially confined to a central region of the domain with no siderophore or siderophore–iron complexes. Thus, the initial data is $$\begin{equation} \left[ \rho, \rho^\prime, n, s_i, s_e, I, C, V \right] = \left\{ \begin{array}{ll} \left[ \rho_0, 0, n_0, s_{i_0}, s_{e_0}, I_0, 0, 0 \right] & \mbox{if}\ r<R_{plug}, \\ & \\ \left[0, 0, 0, 0, s_{e_{0}}, I_0, 0, 0 \right] & {\mbox{otherwise}}, \end{array} \right. \end{equation}$$(2.2) where |$r$| denotes the distance from the centre of the domain (i.e. the inoculation site) while zero-flux boundary conditions are applied on the boundary |$r=R_{dish}$| for all model variables. 2.2 Numerical solutions The model equations (2.1) with initial data (2.2) were solved using Comsol Multiphysics. Parameter values and initial data were used from the calibrations in Boswell et al. (2002, 2003); Perez-Meranda et al. (2007); Eberi & Collinson (2009); Leventhal et al. (2019) while reasonable assumed values were taken for the complex uptake rate |$r_3$| (Table 2). A typical solution is shown in Fig. 1. Fig. 1. Open in new tabDownload slide Numerical solution of equation (2.1) with initial data (2.2) at times |$t=0.5, 1, 1.5$| (representing days) over a circular domain with representative diameter 4 cm. Parameter values are given in Table 2. For each variable, the colour range is shown as a proportion of their maximum value between times |$t=0$| and |$t=3$| (i.e. when the biomass had collided with the edge of the domain). Fig. 1. Open in new tabDownload slide Numerical solution of equation (2.1) with initial data (2.2) at times |$t=0.5, 1, 1.5$| (representing days) over a circular domain with representative diameter 4 cm. Parameter values are given in Table 2. For each variable, the colour range is shown as a proportion of their maximum value between times |$t=0$| and |$t=3$| (i.e. when the biomass had collided with the edge of the domain). Fig. 2. Open in new tabDownload slide The distance |$r$| at times |$t$| from the centre of the domain where |$\rho +\rho ^\prime$| (dotted), |$C$| (dot-dashed), |$I$| (dashed) and |$V$| (solid) take critical values for differing concentrations of iron. (a) |$I_0=0.0004$|⁠, (b) |$I_0=0.004$|⁠, (c) |$I_0=0.04$| representing one tenth, 1 and 10 times the calibrated initial iron value. (Notice that the complex |$V$| expands as a ‘ring’ and hence its inner and outer extents are shown.) The critical concentrations are defined to be 0.0181 for biomass, |$4\times 10^{-5}$| for iron, 0.0679 for siderophores and |$1.204\times 10^{-4}$| for the complex, representing one tenth of their maximum values for the numerical solution with |$I_0=0.0004$|⁠. Fig. 2. Open in new tabDownload slide The distance |$r$| at times |$t$| from the centre of the domain where |$\rho +\rho ^\prime$| (dotted), |$C$| (dot-dashed), |$I$| (dashed) and |$V$| (solid) take critical values for differing concentrations of iron. (a) |$I_0=0.0004$|⁠, (b) |$I_0=0.004$|⁠, (c) |$I_0=0.04$| representing one tenth, 1 and 10 times the calibrated initial iron value. (Notice that the complex |$V$| expands as a ‘ring’ and hence its inner and outer extents are shown.) The critical concentrations are defined to be 0.0181 for biomass, |$4\times 10^{-5}$| for iron, 0.0679 for siderophores and |$1.204\times 10^{-4}$| for the complex, representing one tenth of their maximum values for the numerical solution with |$I_0=0.0004$|⁠. The model biomass expanded outwards in a radially symmetric manner and was preceded by an increased density of model tips. The external substrate was depleted in regions occupied by the biomass. Indeed, the numerical solutions for these variables mirrored those in Boswell et al. (2003), indicating the validity of the numerical integration scheme utilized in the current study and therefore are consistent both qualitatively and quantitatively with experimental data on the growth of a mycelium in initially uniform nutrient settings (Boswell et al., 2003 and references therein). The siderophore population was greatest at the centre of the domain and expanded beyond the extent of the biomass. The iron distribution was depleted from the middle of the domain outwards and the extent of the depletion exceeded the range of the model biomass. The resultant siderophore–iron complex distribution was greatest in the zone between the model biomass and where the iron population was at its greatest and thus consistent with the complex’s formation where siderophores first encounter iron and where the complex is absorbed by the biomass. While the quantitative concentrations of the siderophores and the siderophore–iron complex cannot be related to experimental data, the depletion of the iron population has the same qualitative features as observed in numerous experiments (e.g. Milagres et al., 1999; Bertrand et al., 2010; Srivastava et al., 2013), namely iron is depleted in a radially symmetric fashion and this depletion extends beyond the extremes of the fungal biomass. Indeed, the formation of the siderophore–iron complex coincides with the depletion of iron, and hence, the extent of the complex |$V$| from the biomass periphery yields information on the magnitude of the zone within which the siderophores operate. Fig. 3. Open in new tabDownload slide The distance |$r$| at times |$t$| from the centre of the domain where |$\rho +\rho ^\prime$| (dotted), |$C$| (dot-dashed), |$I$| (dashed) and |$V$| (solid) take critical values for differing concentrations of external substrate. (a) |$s_{e_0}=0.06$|⁠, (b) |$s_{e_0}=0.3$|⁠, (c) |$s_{e_0}=0.6$| representing one tenth, one half and one multiple of the calibrated value. (Notice that the complex |$V$| expands as a ‘ring’ and hence its inner and outer extents are shown.) The critical concentrations are defined to be 0.01 for biomass, |$4\times 10^{-4}$| for iron, 0.04 for siderophore and |$3.341\times 10^{-4}$| for the complex, representing one tenth of their maximum values for the numerical solution with |$s_{e_0}=0.06$|⁠. Fig. 3. Open in new tabDownload slide The distance |$r$| at times |$t$| from the centre of the domain where |$\rho +\rho ^\prime$| (dotted), |$C$| (dot-dashed), |$I$| (dashed) and |$V$| (solid) take critical values for differing concentrations of external substrate. (a) |$s_{e_0}=0.06$|⁠, (b) |$s_{e_0}=0.3$|⁠, (c) |$s_{e_0}=0.6$| representing one tenth, one half and one multiple of the calibrated value. (Notice that the complex |$V$| expands as a ‘ring’ and hence its inner and outer extents are shown.) The critical concentrations are defined to be 0.01 for biomass, |$4\times 10^{-4}$| for iron, 0.04 for siderophore and |$3.341\times 10^{-4}$| for the complex, representing one tenth of their maximum values for the numerical solution with |$s_{e_0}=0.06$|⁠. 2.2.1 Variations in initial iron concentration The extent of the biomass, siderophores, iron and siderophore–iron complexes depend on the concentration of iron as shown in Fig. 2. The extent is defined to be the boundary where each concentration is equal to a critical level (stated in the figure legend) and since the siderophore–iron complex advances as a ‘ring’ formation, both the inner and outer boundaries of that structure are shown. In all cases, the extent of the biomass increases approximately linearly over time, indicative of a constant growth rate, and this also marginally increases with |$I_0$|⁠, consistent with the use of that resource to further promote growth, and therefore has similar characteristics to other modelling investigations and experimental results (e.g. Prosser & Trinci, 1979, where tip vesicles are analogous to internal substrate). The siderophore and iron extent both decline with increasing iron concentration because of the concomitant increased rate of complex formation. Consequently, the extent of the complexes increases with the initial iron concentration |$I_0$| so that whereas for reduced initial iron concentrations, the complexes are only found in the vicinity of the biomass edge (Fig. 2(a)), for greater concentrations the complexes are found throughout most of the domain (e.g. Fig 2(c)). This observation significantly extends experimental results that focus only on the uptake of iron from the growth medium and, due to the difficulties in tracking siderophores, not their distribution prior to or after forming the iron complexes. Throughout Fig. 2, the sudden increase in the extent of the siderophore and siderophore–iron complex populations close to |$r=2$| is due to their interactions with the boundary at |$r=2$|⁠. 2.2.2 Variations in initial external substrate concentration The extent of the biomass, siderophores, iron and siderophore–iron complex are strongly influenced by the concentration of the external substrate as shown in Fig. 3. Firstly, the extent of the biomass increased with the external substrate due to the increased uptake, branching and model tip extension associated with that resource with the least external substrate corresponding to minimal biomass expansion (Fig. 3(a)), consistent with widely reported data relating fungal growth and productivity to nutrient availability (e.g. Suberkropp, 2011). Additionally, since siderophores are produced at a rate proportional to the internal substrate, the siderophore extent also increased with the external substrate, with the initial internal substrate concentration responsible for an initial but not sustained production of siderophores under reduced external substrate concentrations (Fig. 3(a)). The depletion of the iron increased with the external substrate but not linearly; a ten-fold reduction from the default value of |$s_{e_{0}}$| (Fig. 3(c)) did not result in a ten-fold reduction of the extent of iron (Fig. 3(a)). However, the distribution of the siderophore–iron complexes displayed a highly irregular association with the external substrate. For low concentrations of the external substrate, the complex distribution arose as a narrow ‘ring’, a significant distance away from the biomass periphery (Fig. 3(a)). However, as the external substrate increased, the width of this ‘ring’ increased through a reduction in its inner radius (Fig. 3(b)). As the external substrate concentration increased still further, the inner radius of the ‘ring’ expanded and thus reduced the region of the domain in which the complexes were greatest in concentration (Fig. 3(c)). This nonlinear change in siderophore–iron complex distributions due to external substrate concentrations is likely because of associated variations in the production of siderophores coupled with the formation of the complexes and their subsequent uptake by the biomass. For large concentrations of the external substrate, not only were large amounts of siderophores produced, but also the biomass expanded quickly that enabled a more rapid uptake of siderophore–iron complexes. On the other hand, for reduced concentrations of the external substrate, fewer siderophores were produced, the production of siderophore–iron complexes was thus reduced and their subsequent uptake by the biomass was delayed since biomass expansion was slower. 2.2.3 Cumulative iron uptake As previously explained, microorganisms produce siderophores to acquire iron only when in an iron-deficient state. Consequently, quantitative predictions on the amount of iron obtained by the biomass through the acquisition of the iron–siderophore complexes is fundamental in this model. (Indeed, when the internalized iron concentration reaches such a critical level then siderophore production is ceased.) It has previously been shown that the extent of biomass, siderophores and siderophore–iron complexes depends on the initial concentration of iron and external substrate that will therefore also impact on the ultimate uptake of iron by the biomass. Due to the structure of the model equations, iron is either free in the external environment, combined as complexes with siderophores or has been taken up by the model biomass. Consequently, the cumulative amount of iron acquired by the model biomass at time |$t$| can be easily calculated by considering the difference between the initial iron population and the amount of iron at time |$t$| existing in either their free form (i.e. denoted by |$I$|⁠) or that currently held in complexes (i.e. denoted by |$V$|⁠): $$\begin{equation} \mbox{cumulative iron uptake by time}\ t = \int_{\varOmega} I_0(x,y) \, \mbox{d} \varOmega - \int_{\varOmega} I(x,y,t) \, \mbox{d} \varOmega - \int_{\varOmega} V(x,y,t) \, \mbox{d} \varOmega, \end{equation}$$(2.3) where |$\varOmega$| denotes the entire domain (i.e. the region inside the Petri dish). The cumulative amount of iron obtained by the biomass depended upon the initial amount of iron in the external environment and on the external substrate (Fig. 4). In all instances, there was a sudden increase in the quantity of internally held iron and the rate of increase subsequently declined until boundary effects impacted on this process (approximately at time |$t=2$| for the simulations with large values of |$I_0$| and |$s_{e_0}$|⁠). While there appears to be a near linear relationship between the amount of iron in the external environment and that subsequently obtained by the biomass (Fig. 4(a)), there is a more complex non-linearity between the external substrate concentrations and the amount of iron obtained where a ten-fold reduction in external resources only approximately halves the total amount of iron acquired by the biomass (Fig. 4(b)). Fig. 4. Open in new tabDownload slide The cumulative amount of iron obtained by the biomass. (a) The initial iron concentration is varied where |$I_0=0.004$| (solid line), |$I_0=0.04$| (dashed line), |$I_0=0.0004$| (dotted line) and in all cases |$s_{e_0}=0.6$|⁠. (b) The initial external substrate is varied where |$s_{e_0}=0.6$| (solid line), |$s_{e_0}=0.3$| (dashed line), |$s_{e_0}=0.06$| (dotted line) and in all cases |$I_0=0.004$|⁠. Fig. 4. Open in new tabDownload slide The cumulative amount of iron obtained by the biomass. (a) The initial iron concentration is varied where |$I_0=0.004$| (solid line), |$I_0=0.04$| (dashed line), |$I_0=0.0004$| (dotted line) and in all cases |$s_{e_0}=0.6$|⁠. (b) The initial external substrate is varied where |$s_{e_0}=0.6$| (solid line), |$s_{e_0}=0.3$| (dashed line), |$s_{e_0}=0.06$| (dotted line) and in all cases |$I_0=0.004$|⁠. 3. Siderophore-complex distributions: an algebraic approach The analysis in the previous section essentially focussed on the temporal change in the distances over which the siderophores operated and generated good qualitative agreement with experimental observations (Milagres et al., 1999; Srivastava et al., 2013). In polluted terrestrial environments, combinations of heavy metals and other toxins may be present that inhibit the growth of a fungus (Fomina et al., 2000), in addition to toxicity from the HDTMA visual indicator. In such cases, while siderophores are still released by fungi in an iron-depleted state, the mycelium does not necessarily expand due to the presence of pollutants. Since the standard experimental approach to observing siderophore dynamics relates to observing the reduction in iron from the growth medium, if the initial iron concentration is sufficiently high then small losses may not be visually observable. Here, a simplification of model equation (2.1) is used to construct quantitative predictions on siderophore and siderophore–iron complexes in such settings. It is assumed that a circular biomass in an iron-depleted state is positioned inside a toxic region that prohibits its subsequent expansion. This could represent a situation where a fungus is introduced to a domain exhibiting large concentrations of heavy metals, which, for example, arises in bio-remediation applications. Distributions of iron are positioned outside of the toxic region and therefore siderophores provide the sole means of the biomass obtaining iron. See Fig. 5 for a schematic illustration. A key aspect of this investigation is the distance between the biomass, where the siderophores are produced, and the iron resource, where the complexes are formed. Thus, the radius of the biomass, |$R_1$|⁠, and the distance of the iron from the centre of the biomass, |$R_2$|⁠, are crucial parameters. Fig. 5. Open in new tabDownload slide The circular fungal biomass has radius |$R_1$| centred at the origin while the iron is contained outside the circular region of radius |$R_2$|⁠. The siderophores and the siderophore–iron complexes exist in the ‘ring’ between the two circles. Siderophores are released from the biomass at |$r=R_1$| are converted into complexes at |$r=R_2$|⁠, and the complexes are subsequently taken up by the biomass at |$r=R_1$|⁠. Fig. 5. Open in new tabDownload slide The circular fungal biomass has radius |$R_1$| centred at the origin while the iron is contained outside the circular region of radius |$R_2$|⁠. The siderophores and the siderophore–iron complexes exist in the ‘ring’ between the two circles. Siderophores are released from the biomass at |$r=R_1$| are converted into complexes at |$r=R_2$|⁠, and the complexes are subsequently taken up by the biomass at |$r=R_1$|⁠. The biomass is assumed to release siderophores at a constant rate that subsequently diffuse. Since it is reasonable to assume the diffusive time scale is greater than the reactive time scale, once the siderophores encounter the distribution of iron a siderophore–iron complex is immediately formed and diffuses. When the complex reaches the biomass, it is immediately absorbed so that the iron can be utilized by the biomass. Consequently, the above scenario can be represented using polar coordinates and due to radial symmetry (see also the results in Section 2), there is no variation with the angular coordinate. Thus, consistent with the above approaches, the siderophore population is governed by $$\begin{equation} \frac{\partial C}{\partial t} = \frac{D_c}{r} \frac{\partial}{\partial r} \left( r \frac{\partial C}{\partial r}\right), \qquad \mbox{for}\ \ R_1 < r < R_2. \end{equation}$$(3.1) Since siderophores are released at a constant rate by the biomass and immediately form siderophore–iron complexes once the iron distribution is encountered, the corresponding boundary conditions are $$\begin{equation} D_c \frac{\partial C}{\partial r}(R_1,t) = -k, \qquad \qquad C(R_2,t)=0, \end{equation}$$(3.2) where the flux |$k$| corresponds to the rate siderophores enter the region |$R_1<r<R_2$| from the biomass. It is useful to note that boundary condition (3.2) on |$r=R_1$| is an alternative but eventually equivalent condition obtained from the solution of $$\begin{equation} \begin{split} \frac{\partial C}{\partial t} = & \frac{D_c}{r} \frac{\partial}{\partial r} \left(r \frac{\partial C}{\partial r} \right) + \frac{2k}{R_1} H(R_1-r), \qquad \qquad \mbox{for}\ 0<r<R_2,\\ & \frac{\partial C}{\partial r}(0,t)=0, \qquad \qquad C(R_2,t)=0, \end{split} \end{equation}$$(3.3) where |$H$| denotes the standard Heaviside step function and represents the case where siderophores are produced throughout the region |$r<R_1$| at a constant rate so that after a transient time the flux at |$r=R_1$| is a constant |$-\frac{k}{D_c}$|⁠. For convenience, we use boundary condition (3.2) but will later exploit (3.3) in Section 3.1.2. It is assumed that initially, there are no siderophores in the domain, i.e. $$\begin{equation} C(r,0) =0. \end{equation}$$(3.4) The siderophore–iron complex also undergoes diffusion and hence is modelled using $$\begin{equation} \frac{\partial V}{\partial t} = \frac{D_v}{r} \frac{\partial}{\partial r} \left( r \frac{\partial V}{\partial r}\right), \qquad \mbox{for}\ \ R_1<r<R_2. \end{equation}$$(3.5) The siderophore–iron complex forms immediately upon interaction between the siderophores and the iron distribution. Half the resultant complex continues to diffuse in the outward direction while the other half diffuses back towards the biomass whereupon it is immediately absorbed. Thus, the boundary conditions are given by $$\begin{align} V(R_1,t)=0, \qquad \qquad D_v \frac{\partial V}{\partial r}(R_2,t) = -\frac{1}{2}D_c \frac{\partial C}{\partial r}(R_2,t). \end{align}$$(3.6) It is assumed that at time |$t=0$|⁠, there are no siderophore–iron complexes in the domain and so the initial data is $$\begin{equation} V(r,0)=0. \end{equation}$$(3.7) Notice that equations (3.1) and (3.5) are similar to annihilation models (e.g. Ben-Haim & Redner, 1992) except the annihilation arises from a boundary condition rather than a reaction. It is advantageous to nondimensionalize the model equations before constructing their solution. By introducing |$t^* = \frac{D_v}{R_1^2}t, r^{*} = \frac{r}{R_1}, R=\frac{R_2}{R_1}, D = \frac{D_c}{D_v}, C^* = \frac{D_v C}{k R_1}$| and |$V^*=\frac{2D_v V}{k R_1}$|⁠, the model equations reduce to $$\begin{equation} \frac{\partial C}{\partial t} =\frac{D}{r} \frac{\partial}{\partial r} \left( r \frac{\partial C}{\partial r}\right), \qquad \qquad \mbox{for}\ \ 1<r<R, \end{equation}$$(3.8a) $$\begin{equation} \frac{\partial V}{\partial t} =\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial V}{\partial r}\right), \qquad \qquad \mbox{for}\ \ 1<r<R, \end{equation}$$(3.8b) with boundary conditions and initial data $$\begin{equation}\hskip-60pt D \frac{\partial C}{\partial r}(1,t) = - 1, \qquad V(1,t)=0, \end{equation}$$(3.9a) $$\begin{equation} C(R,t)=0, \qquad D \frac{\partial C}{\partial r}(R,t) = -\frac{\partial V}{\partial r}(R,t), \end{equation}$$(3.9b) $$\begin{equation}\hskip-85pt C(r,0)=V(r,0)=0 \end{equation}$$(3.9c) and where |$^{*}s$| have been dropped for notational convenience. The solutions of equation (3.8) with the initial data and boundary conditions (3.9) are in Appendix A shown to be $$\begin{equation} \begin{split} C(r,t) &= \frac{1}{D} \ln \left(\frac{R}{r}\right) + \sum_{n=1}^{\infty} A_n \phi_n(r) e^{-\lambda_n D t}, \\ V(r,t) &= \ln\left(r\right) -D \sum_{n=1}^{\infty} A_n \frac{\phi_n^\prime(R)}{\psi_n^\prime(R)} \psi_n(r) e^{-\lambda_n D t} + \sum_{n=1}^{\infty} E_n \omega_n(r) e^{-\mu_n t}, \end{split} \end{equation}$$(3.10) where the eigenvalues |$\lambda _n$| and |$\mu _n$| are the roots of the characteristic equations $$\begin{equation}\begin{split} \frac{J_1\left(\sqrt{\lambda_n}\right)}{J_0\left(\sqrt{\lambda_n}R\right)} - \frac{Y_1\left(\sqrt{\lambda_n}\right)}{Y_0\left(\sqrt{\lambda_n} R\right)} &=0,\\ \frac{J_1\left(\sqrt{\mu_n} R\right)}{J_0\left(\sqrt{\mu_n}\right)} - \frac{Y_1\left(\sqrt{\mu_n} R\right)}{Y_0\left(\sqrt{\mu_n}\right)} &=0, \end{split} \end{equation}$$(3.11) respectively, the eigenfunctions |$\phi _n(r)$|⁠, |$\psi _n(r)$| and |$\omega _n(r)$| are given by $$\begin{equation*}\begin{split} \phi_n(r) &= \frac{Y_0\left(\sqrt{\lambda_n} r\right)}{Y_0\left(\sqrt{\lambda_n}R\right)} - \frac{J_0\left(\sqrt{\lambda_n}r\right)}{J_0\left(\sqrt{\lambda_n} R\right)},\\ \psi_n(r) &= \frac{Y_0\left(\sqrt{\lambda_n D} r\right)}{Y_0\left(\sqrt{\lambda_n D}\right)} - \frac{J_0\left(\sqrt{\lambda_n D} r\right)}{J_0\left(\sqrt{\lambda_n D}\right)}, \\ \omega_n(r) &= \frac{Y_0\left(\sqrt{\mu_n} r\right)}{Y_0\left(\sqrt{\mu_n}\right)} - \frac{J_0\left(\sqrt{\mu_n} r\right)}{J_0\left(\sqrt{\mu_n}\right)}, \end{split} \end{equation*}$$ where |$J_m$| and |$Y_m$| (⁠|$m=0,1$|⁠) are the Bessel functions of the first and second kind, respectively, and |$A_n$| and |$E_n$| are given by $$\begin{equation*}\begin{split} A_n &=\frac{2 \phi_n(1)}{D \left[\phi_n^2(1)\lambda_n - R^2 (\phi^\prime_n(R))^2\right]}, \\ E_n&= \frac{2 \omega_n(R) \left[ 1 - R \mu_n D \sum_{m=1}^{\infty} \frac{A_m \phi_m^\prime(R)}{\mu_n - D \lambda_m} \right]}{ \left(\omega_n^\prime(1)\right)^2 - \mu_n R^2 \omega_n^2(R)}, \end{split} \end{equation*}$$ provided |$\mu _n \neq D \lambda _m$| for all eigenvalues |$\mu _n$| and |$\lambda _m$|⁠. Notice that the eigenvalues |$\lambda _n$| and |$\mu _n$| from equation (3.11) correspond to the zeros of a cross product of Bessel functions, which have long been studied (e.g. Fettis & Caslin, 1966). 3.1 Results using numerically computed eigenvalues The eigenvalues in equation (3.11) were computed numerically in Matlab enabling the calculation of solutions in (3.10). Since these solutions involve generalized Fourier series, in the investigations below, the summations in equation (3.10) are truncated after 10 terms since the inclusion of further terms produced graphically indistinguishable results. 3.1.1 Typical results The temporal changes in the distributions of the siderophore and siderophore–iron complexes, as obtained from equation (3.10), are shown in Fig. 6. Initially, both distributions are zero throughout the domain. Due to the influx of siderophores at the |$r=1$| boundary, the distribution of siderophores increases accordingly. After a sufficient time has passed, the siderophore population has extended across the domain to reach the |$r=R$| boundary. Accordingly, the production of the siderophore–iron complexes is initiated and this continues to increase so that the complexes subsequently diffuse back across the domain where they are absorbed at the |$r=1$| boundary. As expected, the siderophore distribution approaches its steady state prior to that of the siderophore–iron complex distribution and where the steady states |$C_S(r)$| and |$V_S(r)$| are respectively given by the leading terms in equation (3.10), i.e. $$\begin{equation} C_S(r) = \frac{1}{D} \ln \left(\frac{R}{r}\right), \qquad \qquad V_S(r) = \ln(r), \qquad \qquad \mbox{for}\ 1<r<R. \end{equation}$$(3.12) Fig. 6. Open in new tabDownload slide The distribution of siderophores |$C$| (solid lines) and siderophore–iron complexes |$V$| (dashed lines) from equation (3.10) with |$D=1$| and |$R=2$| are shown at times |$t=0, 0.2, 0.4,\ldots , 2$|⁠. Fig. 6. Open in new tabDownload slide The distribution of siderophores |$C$| (solid lines) and siderophore–iron complexes |$V$| (dashed lines) from equation (3.10) with |$D=1$| and |$R=2$| are shown at times |$t=0, 0.2, 0.4,\ldots , 2$|⁠. Thus increases in |$R$|⁠, representing the relative difference between the location of the iron and the extent of the biomass, results in increases in the density of both the siderophore and the siderophore–iron complex throughout the domain. 3.1.2 Siderophore-complex distribution: numerical predictions The above algebraic solution can be compared to the numerical solutions of model equation (2.1). To simulate the configuration in Fig. 5, |$R_{dish}$| was taken to be 0.45 and the iron was located in the region within a distance 0.05 of the boundary of that boundary so that |$R_2=0.4$| and |$R_1=R_{plug}=0.2$| as before. To represent large concentrations of iron, the calibrated value of |$I_0$| in Table 2 was increased 100 fold and was assumed to be continually replenished upon the production of siderophore–iron complexes and was implemented by removing the corresponding depletion term in equation (2.1f). Finally, the biomass was prevented from expanding from its initial distribution by setting both |$v$| and |$D_n$| to be zero. To best compare the output of the full model equations (2.1) to the algebraic solutions (3.10), notice that the siderophore population in equation (2.1g) can be approximated by (3.3) and that after the same nondimensionalization described above the steady-state solution is $$\begin{equation*} C(r) = \begin{cases}\frac{1+2\ln(R)-{r}^2}{2D},& {\mbox{for}}\ 0<r<1, \\[6pt] \frac{1}{D}\ln\left( \frac{R}{r}\right), & {\mbox{for}}\ 1<r<R, \end{cases} \end{equation*}$$ and therefore satisfies |$D\frac{\partial C}{\partial r}(1)=-1$|⁠, consistent with equation (3.9a). Using this approach, it is seen that the siderophore and siderophore–iron complex populations develop in a similar way to that seen previously (Fig. 7). Indeed, the main difference between the numerical and algebraic solutions arises at |$r=1$| for small times due to the immediate uptake of siderophore–iron complex in the latter (via boundary condition (3.9a)) compared to a more prolonged process in the former (represented by the reaction term in equation (2.1h)). Thus, there is clearly a strong qualitative and quantitative agreement between the algebraic and numerical solutions. Fig. 7. Open in new tabDownload slide The distribution of nondimensionalized siderophores |$C$| (solid lines) and siderophore–iron complexes |$V$| (dashed lines) from equation (2.1) vary over the distance |$r$| from the centre of the domain and the densities increase over time. Except for key parameters described in text, parameter values are given in Table 2 and the distributions are shown at times |$t=0, 0.2, 0.4,\ldots , 2$|⁠. Fig. 7. Open in new tabDownload slide The distribution of nondimensionalized siderophores |$C$| (solid lines) and siderophore–iron complexes |$V$| (dashed lines) from equation (2.1) vary over the distance |$r$| from the centre of the domain and the densities increase over time. Except for key parameters described in text, parameter values are given in Table 2 and the distributions are shown at times |$t=0, 0.2, 0.4,\ldots , 2$|⁠. 3.1.3 Dependence on parameter values Having demonstrated the agreement between the numerical solution of the full set of PDEs (2.1) and the reduced versions (3.8) with boundary conditions (3.9), the algebraic solutions of the reduced equations can be used to obtain useful predictions on the temporal behaviour of the siderophore–iron interactions. While the ultimate effect of the parameters on the final steady-state distribution of the siderophores and the complexes are obvious through equation (3.12), their involvement in the time taken to reach their stationary distributions is less clear. To illustrate the delay in approaching the equilibrium distributions |$C_S(r)$| and |$V_S(r)$|⁠, consider the normalized functions $$\begin{equation} Q_C(t) = \frac{\int_{1}^{R}r C(r,t) \, \mbox{d}r}{\int_{1}^{R} r C_S(r) \, \mbox{d}r}, \qquad \qquad Q_V(t) = \frac{\int_{1}^{R}r V(r,t) \, \mbox{d}r}{\int_{1}^{R} r V_S(r) \, \mbox{d}r}, \end{equation}$$(3.13) which at time |$t=0$| take a value of 0 and approach 1 as the respective distributions approach their equilibria, and therefore represent the ratios of the total amount of each population to their final amount. Notice that both numerators and denominators of (3.13) can be calculated using integration by parts. Figure 8 illustrates the convergence of the siderophore and siderophore–iron complex to their equilibrium distributions for different values of |$D$|⁠. In all cases, the siderophore distribution approaches its equilibria in advance of the siderophore–iron complex. As |$D$| increases, the siderophore distribution approaches its equilibrium more rapidly, consistent with the corresponding increase in movement rates for that population. The delay between the siderophore and complex distributions approaching their equilibria increases with |$D$| up to a limiting value. For |$D\ll 1$|⁠, there is a noticeable lag period before |$Q_V(t)$| increases, corresponding to the time taken for the siderophores to reach the |$r=R$| boundary and initiate the formation of the siderophore–iron complexes; when |$D\gg 1$| no such lag is present due to the comparative reduction in transit time between the two boundaries. Fig. 8. Open in new tabDownload slide The functions |$Q_C(t)$| (solid) and |$Q_V(t)$| (dashed) with |$R=2$| are shown for (a) |$D=0.1$|⁠, (b) |$D=1$|⁠, (c) |$D=10$| and (d) |$D=100$|⁠. Note the different time scale in (a). Fig. 8. Open in new tabDownload slide The functions |$Q_C(t)$| (solid) and |$Q_V(t)$| (dashed) with |$R=2$| are shown for (a) |$D=0.1$|⁠, (b) |$D=1$|⁠, (c) |$D=10$| and (d) |$D=100$|⁠. Note the different time scale in (a). Fig. 9. Open in new tabDownload slide Absolute differences in the times taken by |$Q_v(r,t)$| to approach |$\hat{Q}=0.9$| (solid) and |$\hat{Q}=0.99$| (dashed) obtained using equation (3.10) with numerically computed eigenvalues |$\lambda _n, \mu _n$| for |$n=1,\ldots 10$| and approximation (3.14) using only |$\mu _1$| and |$\lambda _1$| with |$R=2$| for different values of |$D$|⁠. Fig. 9. Open in new tabDownload slide Absolute differences in the times taken by |$Q_v(r,t)$| to approach |$\hat{Q}=0.9$| (solid) and |$\hat{Q}=0.99$| (dashed) obtained using equation (3.10) with numerically computed eigenvalues |$\lambda _n, \mu _n$| for |$n=1,\ldots 10$| and approximation (3.14) using only |$\mu _1$| and |$\lambda _1$| with |$R=2$| for different values of |$D$|⁠. For large values of |$D$|⁠, |$Q_C(t)$| very quickly approaches its equilibrium value of unity so that |$C(r,t)$| can be approximated by its steady-state distribution |$C_S(r)$| while an asymptotic expression can be constructed for the distribution |$V(r,t)$| by taking leading order terms so that $$\begin{equation*} V \approx \ln(r)+E_1 \omega_1(r) e^{-\mu_1 t}. \end{equation*}$$ Consequently for large |$D$|⁠, by noting |$r\omega _1(r) = -\frac{1}{\mu _1} \left (r \omega _1^{\prime }\right )^\prime$| and using integration by parts, it follows that $$\begin{equation*} Q_V(t) \approx 1 - \frac{4 E_1 e^{-\mu_1 t} \left[ R \omega^\prime_1(R) - \omega^\prime_1(1) \right]} {\mu_1 \left[ 1 - R^2 + 2R^2 \ln R \right]} \end{equation*}$$ and therefore the approximate time |$\hat{t}$| for |$Q_V(t)$| to obtain a value |$\hat{Q}$| for large |$D$| is given by $$\begin{equation} \hat{t} = \frac{1}{\mu_1} \ln \left( \frac{ 4 E_1 \left[ R \omega_1^\prime(R) - \omega^\prime_1(1)\right]} { \left( 1- \hat{Q} \right) \mu_1 \left( 1 - R^2 + 2R^2 \ln R \right)} \right). \end{equation}$$(3.14) Notice that the coefficient |$E_1$| involves a summation of terms including |$\lambda _n$|⁠. However, the asymptotic approximation in equation (3.14) where the coefficient |$E_1$| has been truncated to only the leading term involving |$\lambda _1$| agrees well with solutions obtained by algebraically solving equation (3.10) with the first 10 eigenvalues of both |$\lambda _n$| and |$\mu _n$| using Matlab and their differences reduce with increasing |$D$| (Fig. 9). Variations in the domain size |$R$| altered the convergence times of the distributions |$Q_c(t)$| and |$Q_v(t)$| to their equilibrium values (Fig. 10). The convergence times for |$Q_c(t)$| and |$Q_v(t)$|⁠, at least for large |$R$|⁠, can be approximated from the corresponding leading eigenvalues, i.e. are given by |$1/\sqrt{\lambda _1}$| and |$1/\sqrt{\mu _1}$|⁠, respectively. Expansions for |$\lambda _1$| and |$\mu _1$| are detailed below (Section 3.3) and consequently the convergence times for |$Q_c(t)$| and |$Q_v(t)$| scale with |$R$| and |$R\sqrt{\ln R}$|⁠, respectively. Fig. 10. Open in new tabDownload slide |$Q_c(t)$| (solid) and |$Q_v(t)$| (dashed) with |$D=1$| are shown for domain sizes (a) |$R=1.5$|⁠, (b) |$R=2$|⁠, (c) |$R=3$| and (d) |$R=6$|⁠. Fig. 10. Open in new tabDownload slide |$Q_c(t)$| (solid) and |$Q_v(t)$| (dashed) with |$D=1$| are shown for domain sizes (a) |$R=1.5$|⁠, (b) |$R=2$|⁠, (c) |$R=3$| and (d) |$R=6$|⁠. 3.2 Approximations of eigenvalues |$\lambda _n$| and |$\mu _n$| The previous algebraic results required the numerical computation of the eigenvalues |$\lambda _n$| and |$\mu _n$| from equation (3.11). A number of authors have constructed algebraic approximations of various Bessel functions (e.g. Bowman, 2003) but these do not immediately help deduce the roots of (3.11) and hence the required eigenvalues. However, by taking an asymptotic series expansion (see Appendix B), approximations to the eigenvalues can be made resulting in an entirely algebraic solution for equation (3.10) under appropriate limits. Indeed, by defining $$\begin{equation*} \begin{split} P_n(p, q) &= \left| \frac{\pi (n-\frac{1}{2})}{q-p} \right|, \\ Q_1(p,q) &= \frac{p+3q}{8 p q (q-p)}, \\ Q_3(p,q) & = \frac{25 p^4-31 p^3 q-36 p^2 q^2+9 p q^3-63 q^4}{384 (q-p)^2 q^3 p^3}, \\ Q_5(p,q) & = \frac{3219p^7 - 6938p^6q + 2279p^5 q^2 + 2040p^4q^3+360p^3q^4+4797p^2q^5-7614pq^6+5697q^7}{15360 p^5 q^5 (q-p)^3}, \end{split} \end{equation*}$$ and provided |$n \gg R$|⁠, the square roots of the eigenvalues for large values of |$R$| can conveniently be expressed as the series $$\begin{equation} \begin{split} \sqrt{\lambda_n} &= P_n(1,R) + \frac{Q_1(1,R)}{P_n(1,R)} + \frac{Q_3(1,R)}{P_n^3(1,R)} + \frac{Q_5(1,R)}{P_n^5(1,R)} + \cdots \\ \sqrt{\mu_n} &= P_n(R,1) + \frac{Q_1(R,1)}{P_n(R,1)} + \frac{Q_3(R,1)}{P_n^3(R,1)} + \frac{Q_5(R,1)}{P_n^5(R,1)} + \cdots \end{split} \end{equation}$$(3.15) By defining the zeroth order approximation as comprising only the first term in the series, the first order approximation comprising only the first two terms and so on, even second order approximations are in close agreement with numerically computed values for all but the smallest eigenvalues |$\lambda _1$| and |$\mu _1$| for |$R=2$| (Table 3). Indeed, good approximations for the eigenvalues |$\lambda _1$| and |$\mu _1$| arise provided sufficient terms in the series approximation are included. Notice, however, that in order to use these approximations, it is necessary that |$D \neq \mu _n / \lambda _m$| for all |$n,m$| to ensure that |$E_n$| is defined. Thus, for example, the zeroth order approximation cannot be used for |$D=1$| (but the first and higher order approximations can still be used). Table 3 Comparison of numerical and analytical values of eigenvalues with |$R=2$| using the approximations in equation (3.15) of stated order. . |$n$| . Numerical . 0th . 1st . 2nd . 3rd . 1 1.7940 1.5708 1.8493 1.7555 1.8440 2 4.8021 4.7124 4.8052 4.8018 4.8021 |$\sqrt{\lambda _n}$| 3 7.9090 7.8540 7.9097 7.9089 7.9090 4 11.0351 10.9956 11.0354 11.0351 11.0351 |$\vdots$| 1 1.3608 1.5708 1.3719 1.3687 1.3504 2 4.6459 4.7124 4.6461 4.6460 4.6459 |$\sqrt{\mu _n}$| 3 7.8142 7.8540 7.8142 7.8142 7.8142 4 10.9671 10.9956 10.9672 10.9671 10.9671 |$\vdots$| . |$n$| . Numerical . 0th . 1st . 2nd . 3rd . 1 1.7940 1.5708 1.8493 1.7555 1.8440 2 4.8021 4.7124 4.8052 4.8018 4.8021 |$\sqrt{\lambda _n}$| 3 7.9090 7.8540 7.9097 7.9089 7.9090 4 11.0351 10.9956 11.0354 11.0351 11.0351 |$\vdots$| 1 1.3608 1.5708 1.3719 1.3687 1.3504 2 4.6459 4.7124 4.6461 4.6460 4.6459 |$\sqrt{\mu _n}$| 3 7.8142 7.8540 7.8142 7.8142 7.8142 4 10.9671 10.9956 10.9672 10.9671 10.9671 |$\vdots$| Open in new tab Table 3 Comparison of numerical and analytical values of eigenvalues with |$R=2$| using the approximations in equation (3.15) of stated order. . |$n$| . Numerical . 0th . 1st . 2nd . 3rd . 1 1.7940 1.5708 1.8493 1.7555 1.8440 2 4.8021 4.7124 4.8052 4.8018 4.8021 |$\sqrt{\lambda _n}$| 3 7.9090 7.8540 7.9097 7.9089 7.9090 4 11.0351 10.9956 11.0354 11.0351 11.0351 |$\vdots$| 1 1.3608 1.5708 1.3719 1.3687 1.3504 2 4.6459 4.7124 4.6461 4.6460 4.6459 |$\sqrt{\mu _n}$| 3 7.8142 7.8540 7.8142 7.8142 7.8142 4 10.9671 10.9956 10.9672 10.9671 10.9671 |$\vdots$| . |$n$| . Numerical . 0th . 1st . 2nd . 3rd . 1 1.7940 1.5708 1.8493 1.7555 1.8440 2 4.8021 4.7124 4.8052 4.8018 4.8021 |$\sqrt{\lambda _n}$| 3 7.9090 7.8540 7.9097 7.9089 7.9090 4 11.0351 10.9956 11.0354 11.0351 11.0351 |$\vdots$| 1 1.3608 1.5708 1.3719 1.3687 1.3504 2 4.6459 4.7124 4.6461 4.6460 4.6459 |$\sqrt{\mu _n}$| 3 7.8142 7.8540 7.8142 7.8142 7.8142 4 10.9671 10.9956 10.9672 10.9671 10.9671 |$\vdots$| Open in new tab When used in (3.10), approximations (3.15) produce results consistent with the full algebraic solutions and are in strong qualitative agreement for small |$R$| (Fig. 11), especially at larger times. Such a result is unsurprising since the approximations in (3.15) were derived from asymptotic expansions of |$J_\nu (z)$| and |$Y_\nu (z)$| for large |$z$| and hence are most applicable for the calculation of |$\lambda _n$| and |$\mu _n$| for large |$n$| (see also Table 3) but the smallest eigenvalues |$\lambda _1$| and |$\mu _1$| exert the greatest influence on the solutions in equation (3.10), particularly at small times. Fig. 11. Open in new tabDownload slide (a) Solution of equation (3.10) with |$R=2$|⁠, |$D=1$| using the eigenvalues computed from equation (3.11), (b), (c) and (d) using the first, second and third approximations from equation (3.15). Profiles of the siderophore distribution |$C$| (solid lines) and the siderophore–iron complex |$V$| (dashed lines) are shown at times |$t=0,0.2,0.4,\ldots 2$|⁠. Fig. 11. Open in new tabDownload slide (a) Solution of equation (3.10) with |$R=2$|⁠, |$D=1$| using the eigenvalues computed from equation (3.11), (b), (c) and (d) using the first, second and third approximations from equation (3.15). Profiles of the siderophore distribution |$C$| (solid lines) and the siderophore–iron complex |$V$| (dashed lines) are shown at times |$t=0,0.2,0.4,\ldots 2$|⁠. 3.3 Approximations of leading eigenvalues |$\lambda _1$| and |$\mu _1$| It was shown above that the approximations (3.15) for |$\lambda _n$| and |$\mu _n$| are least suited for small values of |$n$|⁠, especially |$n=1$|⁠, and also are less suited for large values of |$R$| (see Appendix B). However, the first eigenvalues |$\lambda _1$| and |$\mu _1$| have the most prominent roles in the convergence of the siderophore and siderophore–iron complex to their final steady-state distributions. Hence, an alternative approach to approximating |$\lambda _1$| and |$\mu _1$| is developed here. By observing the behaviour of |$R \sqrt{\lambda _n}$| as |$R \to \infty$|⁠, it follows that (see Appendix C) $$\begin{equation} \sqrt{\lambda_1} = \frac{\zeta_1}{R} - \frac{\pi \zeta_1^2 Y_0(\zeta_1)}{4 R^3 J_{0}^\prime(\zeta_1)} + O\left( R^{-5} \ln(R) \right), \end{equation}$$(3.16) where |$\zeta _1$| is the first root of |$J_0(\zeta )=0$| and is valid for |$\zeta _1 \ll R$|⁠. In a similar way by considering |$R \sqrt{\mu _n}$| as |$R \to \infty$| (Appendix C) $$\begin{equation} \mu_1 = \frac{1}{R^2 \ln (R)} \left[ 2 + \frac{3}{2 \ln (R)} +\frac{5}{6 \ln(R)^2} + O\left(\frac{1}{\ln(R)^3}\right) \right], \end{equation}$$(3.17) which is valid provided |$e \ll R$|⁠. Furthermore, related expressions can be derived for all |$\lambda _n$| and |$\mu _m$| (Appendix C). Table 4 demonstrates that the eigenvalues |$\lambda _1$| and |$\mu _1$| obtained using approximations (3.16) and (3.17) converge to those obtained using the full characteristic equations (3.11) as |$R$| increases. Table 4 Solutions using numerically computed eigenvalues in equation (3.11) are compared to truncating the series to terms in only |$\lambda _1$| and |$\mu _1$| from equations (3.16) and (3.17) for different domain sizes with |$D=1$|⁠. For equations (3.11), |$\bar{t}$| denotes the time taken for |$Q_V(\bar{t})=0.9$|⁠, while the approximation |$\hat{t}$| is obtained from equation (3.14) using eigenvalues (3.16, 3.17). |$t_{C}$| and |$t_{Cf}$| denote the approximate times for the siderophore concentration at |$r=1+0.9(R-1)$| and the flux at |$r=R$| to reach half the steady-state values, respectively (see text for details). . |$R$| . |$\lambda _1$| . |$\mu _1$| . |$\bar{t}$| or |$\hat{t}$| . |$t_{C}$| . |$t_{Cf}$| . Full solution 2 3.2185 1.8517 1.7783 0.3101 0.3118 from (3.11) 6 0.1768 0.0476 57.4074 6.1772 6.2190 20 0.0146 0.0022 1170.4852 78.0076 78.6103 Approx. solutions 2 3.0979 2.1274 1.6872 0.3225 0.3242 using |$\lambda _1$| and |$\mu _1$| 6 0.1776 0.0480 57.1043 6.1651 6.2066 from (3.16) and (3.17) 20 0.0146 0.0022 1177.7715 78.3036 78.9064 . |$R$| . |$\lambda _1$| . |$\mu _1$| . |$\bar{t}$| or |$\hat{t}$| . |$t_{C}$| . |$t_{Cf}$| . Full solution 2 3.2185 1.8517 1.7783 0.3101 0.3118 from (3.11) 6 0.1768 0.0476 57.4074 6.1772 6.2190 20 0.0146 0.0022 1170.4852 78.0076 78.6103 Approx. solutions 2 3.0979 2.1274 1.6872 0.3225 0.3242 using |$\lambda _1$| and |$\mu _1$| 6 0.1776 0.0480 57.1043 6.1651 6.2066 from (3.16) and (3.17) 20 0.0146 0.0022 1177.7715 78.3036 78.9064 Open in new tab Table 4 Solutions using numerically computed eigenvalues in equation (3.11) are compared to truncating the series to terms in only |$\lambda _1$| and |$\mu _1$| from equations (3.16) and (3.17) for different domain sizes with |$D=1$|⁠. For equations (3.11), |$\bar{t}$| denotes the time taken for |$Q_V(\bar{t})=0.9$|⁠, while the approximation |$\hat{t}$| is obtained from equation (3.14) using eigenvalues (3.16, 3.17). |$t_{C}$| and |$t_{Cf}$| denote the approximate times for the siderophore concentration at |$r=1+0.9(R-1)$| and the flux at |$r=R$| to reach half the steady-state values, respectively (see text for details). . |$R$| . |$\lambda _1$| . |$\mu _1$| . |$\bar{t}$| or |$\hat{t}$| . |$t_{C}$| . |$t_{Cf}$| . Full solution 2 3.2185 1.8517 1.7783 0.3101 0.3118 from (3.11) 6 0.1768 0.0476 57.4074 6.1772 6.2190 20 0.0146 0.0022 1170.4852 78.0076 78.6103 Approx. solutions 2 3.0979 2.1274 1.6872 0.3225 0.3242 using |$\lambda _1$| and |$\mu _1$| 6 0.1776 0.0480 57.1043 6.1651 6.2066 from (3.16) and (3.17) 20 0.0146 0.0022 1177.7715 78.3036 78.9064 . |$R$| . |$\lambda _1$| . |$\mu _1$| . |$\bar{t}$| or |$\hat{t}$| . |$t_{C}$| . |$t_{Cf}$| . Full solution 2 3.2185 1.8517 1.7783 0.3101 0.3118 from (3.11) 6 0.1768 0.0476 57.4074 6.1772 6.2190 20 0.0146 0.0022 1170.4852 78.0076 78.6103 Approx. solutions 2 3.0979 2.1274 1.6872 0.3225 0.3242 using |$\lambda _1$| and |$\mu _1$| 6 0.1776 0.0480 57.1043 6.1651 6.2066 from (3.16) and (3.17) 20 0.0146 0.0022 1177.7715 78.3036 78.9064 Open in new tab These approximations for the leading eigenvalues can also be used with the normalized function |$Q_{V}(t)$| in equation (3.13) to estimate the time taken by the siderophore–iron complex to approach its steady-state distribution. Table 4 compares the numerically computed time |$\bar{t}$| such that |$Q_V(\bar{t})=0.9$| using the eigenvalues from the solutions of equation (3.11) to the approximation in equation (3.14) for |$\hat{t}$| where the summation used in |$E_1$| is restricted to its leading term, i.e. that involving only |$\lambda _1$|⁠. The approximations of |$\lambda _1$| and |$\mu _1$| obtained above also allow the derivation of simple expressions relating to the spread of siderophores and the resultant uptake of iron at the |$r=R$| boundary. In particular, by truncating the series to terms only involving |$\lambda _1$| and |$\mu _1$| in the solutions for |$C(r,t)$| in equation (3.10), the approximate time taken |$t_{C}$| for the siderophore density to reach a concentration |$C^\dagger$| at |$r=r^\dagger$| (where |$1<r^\dagger <R$| and |$0<C^\dagger <C_S(r^\dagger )$|⁠) can be shown to satisfy $$\begin{align} t_{C}= - \frac{1}{\lambda_1 D} \ln \left( \frac{DC^\dagger - \ln(R/r^\dagger)}{A_1 D \phi_1(r^\dagger)} \right) \end{align}$$(3.18) while the flux of the siderophores at the boundary |$r=R$| corresponds to the acquisition of iron by the siderophores and the approximate time taken |$t_{Cf}$| for this rate to reach a value |$C_r^\dagger$| (where |$0>C_r^\dagger > -(R D)^{-1}$|⁠, representing the value at equilibrium) satisfies $$\begin{equation} t_{Cf} = -\frac{1}{\lambda_1 D} \ln \left( \frac{ R D C_r^\dagger +1}{ R D A_1 \phi_1^\prime(R)} \right). \end{equation}$$(3.19) These expressions clearly illustrate the effect of the diffusion coefficient |$D$| and the radius |$R$| on the delay until the iron begins to be acquired by the siderophores. In Table 4, the approximations in equations (3.18) and (3.19) using the approximated eigenvalues (3.16) and (3.17) are compared to the corresponding algebraic solutions from equation (3.10) with numerically computed eigenvalues from equation (3.11). The simple approximations using (3.16) and (3.17) are in strong qualitative and quantitative agreement with the full algebraic solution and the agreement improves as |$R$| is increased due to two independent reasons; firstly, the approximation of the leading eigenvalues improves as |$R \to \infty$| and, secondly, as |$R$| increases, it takes longer for the siderophores to reach the exterior boundary at |$r=R$| and hence the second and higher eigenvalues play less significant roles in determining the distributions of |$C(r,t)$| and |$V(r,t)$|⁠. 4. Discussion Siderophores play a central role in how microorganisms acquire important elements. While there are known to be hundreds of different types of siderophores with various functionalities, the most studied relationship is that with iron and thus the subject of this investigation. Indeed, it has recently been shown that siderophores significantly increase the rate at which bacteria acquire this important resource compared to alternative methods (Niehus et al., 2017; Leventhal et al., 2019). Equation (2.1) represents, to the authors’ knowledge, the first mathematical model of iron uptake in fungi mediated through siderophores. The numerical simulations of the model equations display the same qualitative features observed in experiments regarding the extraction of iron from a solid growth medium; specifically, there is a radially symmetric depletion of the iron that extends beyond the edge of the expanding biomass (Fig. 1) and that this region expands initially in an approximately linear fashion at rates determined by local conditions (Figs 2 and 3). In limiting conditions, e.g. Fig. 3(a), the expansion of the siderophore distribution and the concomitant depletion of the iron concentration was clearly less than linear and instead the extent of the iron depletion appeared to increase with the square root of time, consistent with the reduced production and diffusive movement properties of the siderophores. A key feature of the model was its ability to predict the cumulative amount of iron taken up by the biomass through the absorption of the iron–siderophore complexes, as represented by equation (2.3). Such time-dependent data is difficult to obtain experimentally through either direct or indirect means as destructive sampling of the biomass provides the most accurate measurements of the former, while the latter is limited since there is currently no convenient procedure to measure siderophore populations given their diversity. Nonetheless, our model clearly has the potential to make such quantitative predictions on iron acquisition by mycelial fungi. Moreover, further refinements should account for such siderophore diversity and the different pathways through which iron is utilized by fungi following its acquisition (e.g. Howard, 1999). It should also be noted that the model equations represent a simplification of how a combination of different nutrients can impact on the growth and function of a fungal mycelium through the merger of internalized iron and the generic substrate. While alternative approaches have been used to model how fungi utilize combinations of nutrients and essential elements (e.g. Lamour et al., 2000), due to the generalized treatment of the iron pathway once that substance was internalized by the fungus, the precise role of iron on key morphological processes was not isolated in this current study and therefore remains an important avenue for future investigations that would necessitate the inclusion of feedback processes by restricting siderophore production to prevent excessive accumulation of iron. Key features of the numerical solution of the full set of equations (2.1) were captured in the algebraic solutions of the reduced set of equations (3.8), including the constant uptake rate of iron for all but small times. Indeed, there was strong qualitative and quantitative agreement between the full numerical solutions and the algebraic simplifications in the distributions of siderophores and siderophore–iron complexes (Figs 6 and 7). The nondimensionalization used to construct the algebraic solutions (3.10) introduced the parameter |$D$| representing the ratio of the diffusion coefficients of the siderophores and the siderophore–iron complexes. Since the diffusion coefficient of the complexes is less than that of the siderophores (due to obvious differences in their molecular weight), it follows in application that |$D>1$| and therefore siderophores are released and complexes are formed more rapidly than they are acquired by the biomass until equilibrium is reached (Fig. 8). Consequently, equation (3.14) with |$\hat{Q}=0.9$| (or 0.99) is expected to provide a reasonable estimate for the time taken for the siderophore–iron complex distribution to approach its equilibria. The same algebraic solutions also demonstrated the impact of domain size on siderophore and siderophore–iron complex distribution. Specifically, greater distances between the biomass and the source of iron resulted in greater concentrations of both populations (equation (3.12)). An important consequence of the model equations is the ability to calculate the cumulative amount of iron taken up by the biomass through the release of siderophores and the subsequent acquisition of the siderophore–iron complexes. Other than during an initial transient period, the total uptake rate of iron was approximately linear (Fig. 4) except when influenced by boundary effects. Indeed, this same qualitative feature is captured in the reduced model in Section 3 by observing that for large |$D$| (i.e. when |$D_c \gg D_v$|⁠), the uptake of iron corresponded to the flux of the complex at |$r=1$| which to leading order from equation (3.10) is given by $$\begin{equation*} \left. \frac{\partial V}{\partial r} \right|{}_{r=1} \approx 1 + E_1 \omega_1^\prime(1) e^{-\mu_1 t} \end{equation*}$$ and tends to the constant unity. However, this rate was heavily influenced by local conditions. While an increased concentration of external substrate resulted in an increase of iron extracted from the growth domain and internalized by the biomass, the relationship was highly nonlinear; a ten-fold increase in external substrate only approximately doubled the amount of iron obtained by the biomass. However, the observation that external resources can influence the depletion of iron from the growth environment clearly has important consequences in the bio-technological applications of fungi. While the algebraic results presented in this paper have focussed on radial geometry, similar treatments are possible in other domains including a single-dimension Cartesian and spherical radial geometries. (Indeed, by introducing |$x$| so that |$r=R_1+(R_2-R_1) x$| and letting |$R_2-R_1 \to \infty$|⁠, equations (3.1) and (3.5) can be easily converted into a one-dimensional Cartesian geometry with spatial coordinate |$x$| resulting in Fourier series solutions for the siderophore and complex populations. Such a situation has been thoroughly explored in Choudhury (2019).) In our calculations, the algebraic solutions (3.10) are defined provided the nondimensionalized diffusion coefficient |$D$| is not a ratio of the eigenvalues |$\lambda _n$| and |$\mu _m$| for all |$n,m$|⁠. In one-dimensional Cartesian geometry, the equivalent restriction corresponds to |$D$| not being a ratio of squares of odd numbers (however, alternative solutions can be constructed by selecting an alternative form for |$\hat{V}$| in Appendix A, equation (A15)). Moreover, similar issues arise in the spherical radial geometry case. We cannot provide any physical reasoning behind this limitation. Further interesting analysis would concern the implementation of moving boundary conditions consistent with the depletion of the iron concentration and the advancement of the fungal biomass. Such a situation would more closely represent the scenarios considered in Section 2. Siderophores are extensively used by microorganisms to obtain essential metals, in particular iron. In this work, we have constructed and investigated the first mathematical model of their use by fungi. The qualitative behaviour of the model is consistent with known experiments and quantitative predictions have been made on how local conditions influence the amount of iron obtained by the fungus along with how key distributions involving siderophore function change over time. It remains to develop a suitable experimental technique to verify these predictions. We note that the model does not consider how the fungus subsequently uses the iron it has obtained and this is therefore an important challenge for future modelling investigations. References Ahmed , E. & Holmström , S. J. M. ( 2014 ) Siderophores in environmental research: roles and applications . Microb. Biotechnol. , 73 , 196 – 208 . Google Scholar Crossref Search ADS WorldCat Andrews , M. Y. , Santelli , C. M. & Duckworth , O. W. ( 2016a ) Digital image quantification of siderophores on agar plates . Data Brief , 6 , 890 – 898 . Google Scholar Crossref Search ADS WorldCat Andrews , M. Y. , Santelli , C. M. & Duckworth , O. W. ( 2016b ) Layer plate CAS assay for the quantification of siderophore production and determination of exudation patterns for fungi . J. Microbiol. Methods , 121 , 41 – 43 . Google Scholar Crossref Search ADS WorldCat Beard , J. L. ( 2008 ) Why iron deficiency is important in infant development . J. Nutr. , 138 , 2534 – 2536 . Google Scholar Crossref Search ADS PubMed WorldCat Bellenger , J. P. , Wichard , T., Kustka , A. B. & Kraepiel , A. M. L. ( 2013 ) Uptake of molybdenum and vanadium by a nitrogen-fixing soil bacterium using siderophores . Nat. Geosci. , 1 , 243 – 246 . Google Scholar Crossref Search ADS WorldCat Ben-Haim , E. & Redner , S. ( 1992 ) Inhomogeneous two-species annihilation in the steady state . J. Phys. Math. Gen. , 25 , 575 – 583 . Google Scholar Crossref Search ADS WorldCat Bertrand , S. , Bouchara , J. P., Venier , M. C., Richomme , P., Duval , O. & G. Larcher. ( 2010 ) N|$^{\alpha }$|-methyl coprogen B, a potential marker of the airway colonization by Scedosporium apiospermum in patients with cystic fibrosis . Med. Mycol. , 48 , S98 – 107 . Google Scholar Crossref Search ADS PubMed WorldCat Bogumił , A. , Paszt , L. S., Lisek , A., Trzciński , P. & Harbuzov , A. ( 2013 ) Identification of new Trichoderma strains with antagonistic activity against Botrytis cinerea . Folia Hortic. , 25 , 123 – 132 . Google Scholar Crossref Search ADS WorldCat Boswell , G. P. , Jacobs , H., Davidson , F. A., Gadd , G. M. & Ritz , K. ( 2002 ) Functional consequences of nutrient translocation in mycelial fungi . J. Theor. Biol. , 217 , 459 – 477 . Google Scholar Crossref Search ADS PubMed WorldCat Boswell , G. P. , Jacobs , H., Davidson , F. A., Gadd , G. M. & Ritz , K. ( 2003 ) Growth and function of fungal mycelia in heterogeneous environments . Bull. Math. Biol. , 65 , 447 – 477 . Google Scholar Crossref Search ADS PubMed WorldCat Boswell , G. P. , Jacobs , H., Ritz , K., Gadd , G. M. & Davidson , F. A. ( 2007 ) The development of fungal networks in complex environments . Bull. Math. Biol. , 69 , 605 – 634 . Google Scholar Crossref Search ADS PubMed WorldCat Boukhalfa , H. , Lack , J., Reilly , S. D., Herman , L. & Neu , M. P. ( 2003 ) Siderophore production and facilitated uptake of iron and plutonium in P. putida . AIP Conference Proceedings , vol. 673 . New York: American Institute of Physics, pp. 343 – 344 . Google Scholar Crossref Search ADS WorldCat Bowman , F. ( 2003 ) Introduction to Bessel Functions . Dover . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Braud , A. , Hoegy , F., Jezequel , K., Lebeau , T. & Schalk , I. J. ( 2009 ) New insights into the metal specificity of the Pseudomonas aeruginosa pyoverdine-iron uptake pathway . Environ. Microbiol. , 11 , 1079 – 1091 . Google Scholar Crossref Search ADS PubMed WorldCat Choudhury , M. J. A. ( 2019 ) Mathematical modelling of fungal interactions . Ph.D. Thesis , University of South Wales . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Choudhury , M. J. A. , Trevelyan , P. M. J. & Boswell , G. P. ( 2018 ) A mathematical model of nutrient influence on fungal competition . J. Theor. Biol. , 438 , 9 – 20 . Google Scholar Crossref Search ADS PubMed WorldCat Eberi , H. J. & Collinson , S. ( 2009 ) A modelling and simulation study of siderophore mediated antagonism in dual-species biofilms . Theor. Biol. Med. Model. , 6 , 1 – 16 . Google Scholar Crossref Search ADS PubMed WorldCat Fettis , H. E. & Caslin , J. C. ( 1966 ) An extended table of zeros of cross products of Bessel functions . Report No. ARL 66-0023, Aerospace Research Laboratories, Office of Aerospace Research, United States Air Force, Wright-Patterson Air Force Base, Ohio . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Fomina , M. , Ritz , K. & Gadd , G. M. ( 2000 ) Negative fungal chemotropism to toxic metals . FEMS Microbiol. Lett. , 193 , 207 – 211 . Google Scholar Crossref Search ADS PubMed WorldCat Ghosh , S. K. , Pal , S. & Chakraborty , N. ( 2015 ) The qualitative and quantitative assay of siderophore production by some microorganism and effects of different media on its production . Int. J. Chem. Sci. , 13 , 1621 – 1629 . Google Scholar OpenURL Placeholder Text WorldCat Gruhn , C. M. , Gruhn , A. V. & Miller , O. K. ( 1992 ) Boletinellus meruliodes alter root morphology of pinus densiflora without mycorrhizal formation . Mycologia , 84 , 528 – 533 . Google Scholar Crossref Search ADS WorldCat Haas , H. ( 2014 ) Fungal siderophore metabolism with a focus on Aspergillus fumigatus . Nat. Prod. Rep. , 31 , 1266 – 1276 . Google Scholar Crossref Search ADS PubMed WorldCat Harrison , J. ( 2009 ) Fast and accurate Bessel function computation . Proceedings of the 2009 19th IEEE Symposium on Computer Arithmetic . Washington DC : IEEE Computer Society , pp. 104 – 113 . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Howard , D. H. ( 1999 ) Acquisition, transport, and storage of iron by pathogenic fungi . Clin. Micriobiol. Rev. , 12 , 394 – 404 . Google Scholar Crossref Search ADS WorldCat Johnson , L. ( 2008 ) Iron and siderophores in fungal–host interactions . Mycol. Res. , 112 , 170 – 183 . Google Scholar Crossref Search ADS PubMed WorldCat Kraemer , S. M. , Butler , A., Borer , P. & Cervini-Silva , J. ( 2005 ) Siderophores and the dissolution of iron-bearing minerals in marine systems . Rev. Mineral. Geochem. , 59 , 53 – 84 . Google Scholar Crossref Search ADS WorldCat Lamour , A. , van den Bosch , F., Termorshuizen , A. J. & Jeger , M. J. ( 2000 ) Modelling the growth of soil-borne fungi in response to carbon and nitrogen . Math. Med. Biol. , 17 , 329 – 346 . Google Scholar Crossref Search ADS WorldCat Leventhal , G. E. , Ackermann , M. & Schiessl , K. T. ( 2019 ) Why microbes secrete molecules to modify their environment: the case of iron-chelating siderophores . J. R. Soc. Interface , 16 , 20180674 . Google Scholar Crossref Search ADS PubMed WorldCat Machuca , A. & Milagres , A. M. F. ( 2003 ) Use of CAS-agar plate modified to study the effect of different variables on the siderophore production by aspergillus . Lett. Appl. Microbiol. , 36 , 177 – 181 . Google Scholar Crossref Search ADS PubMed WorldCat Marschner , H. ( 1995 ) Mineral Nutrition of Higher Plants . Academic Press . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC McLoughin , T. J. , Quinn , J. P., Bettermann , A. & Bookland , R. ( 1992 ) Pseudomonas cepacia suppression of sunflower wilt fungus and role of antifungal compounds in controlling the disease . Appl. Environ. Microbiol. , 58 , 1760 – 1763 . Google Scholar Crossref Search ADS PubMed WorldCat Milagres , A. M. F. , Machuca , A. & Napoleão , D. ( 1999 ) Detection of siderophore production from several fungi and bacteria by a modification of chrome Azurol S (CAS) agar plate assay . J. Microbiol. Methods , 37 , 1 – 6 . Google Scholar Crossref Search ADS PubMed WorldCat Niehus , R. , Picot , A., Oliveira , N. M., Mitri , S. & Foster , K. R. ( 2017 ) The evolution of siderophore production as a competitive trait . Evolution , 71 , 1443 – 1455 . Google Scholar Crossref Search ADS PubMed WorldCat Perez-Meranda , S. , Cabirol , N., George-Tellez , R., Zamudio-Rivera , L. S. & Fernandez , F. J. ( 2007 ) O-CAS, a fast universal method for siderophore detection . J. Microbiol. Methods , 70 , 127 – 131 . Google Scholar Crossref Search ADS PubMed WorldCat Philpott , C. C. , Leidgens , S. & Frey , A. G. ( 2012 ) Metabolic remodeling in iron-deficient fungi . Biochim. Biophys. Acta , 1823 , 1509 – 1520 . Google Scholar Crossref Search ADS PubMed WorldCat Prosser , J. I. & Trinci , A. P. J. ( 1979 ) A model for hyphal growth and branching . J. Gen. Microbiol. , 111 , 153 – 164 . Google Scholar Crossref Search ADS PubMed WorldCat Renshaw , J. C. , Robson , G. D., Trinci , A. P. J., Wiebe , M. G., Livens , F. R., Collison , D. & Taylor , R. J. ( 2002 ) Fungal siderophores: structures, functions and applications . Mycol. Res. , 106 , 1123 – 1142 . Google Scholar Crossref Search ADS WorldCat Riquelme , M. , Aguirre , J., Bartnicki-García , S., Braus , G. H., Feldbrügge , M, U. Fleig, Hansberg , W., Herrar-Estrella , A., Kämper , J., Kück , U., Mouriño-Pérez , R. R., Takeshita , N. & Fischer , R. ( 2018 ) Fungal morphogenesis, from the polarized growth of hyphae to complex reproduction and infection structures . Microbiol. Mol. Biol. Rev. , 82 , e00068–17 . Google Scholar Crossref Search ADS PubMed WorldCat Saha , M. , Sarkar , S., Sarkar , B., Sharma , B. K., Bhattacharjee , S. & Tribedi , P. ( 2016 ) Microbial siderophores and their potential: a review . Environ. Sci. Pollut. Res. Int. , 23 , 3984 – 2999 . Google Scholar Crossref Search ADS PubMed WorldCat Sasirekha , B. & Srividya , S. ( 2016 ) Siderophore production by Pseudomonas aeruginosa FP6, a biocontrol strain for Rhizoctonia solani and Colletotrichum gloesporioides causing diseases in chilli . Agric. Nat. Resour. , 50 , 250 – 256 . Google Scholar OpenURL Placeholder Text WorldCat Schwyn , B. & Neilands , J. B. ( 1987 ) Universal chemical assay for the detection and determination of siderophores . Anal. Biochem. , 160 , 47 – 56 . Google Scholar Crossref Search ADS PubMed WorldCat Srivastava , M. P. , Tiwari , R. & Sharma , N. ( 2013 ) Effect of different cultural variables on siderophores produced by Trichoderma spp . Int. J. Adv. Res. , 1 , 1 – 6 . Google Scholar OpenURL Placeholder Text WorldCat Suberkropp , K. ( 2011 ) The influence of nutrients on fungal growth, productivity, and sporulation during leaf breakdown in streams . Can. J. Bot. , 73 , 1361 – 1369 . Google Scholar Crossref Search ADS WorldCat Van der Helm , D. & Winkelmann , G. ( 1994 ) Hydroxamates and polycarboxylates as iron transport agents (siderophores) in fungi . Metal Ions in Fungi (G. Winkelmann and D. Winge, eds). Marcel Dekker . Verma , V. C. , Singh , S. K. & Prakash , S. ( 2011 ) Bio-control and plant growth promotion potential of siderophore producing endophytic streptomyces from Azadirachta indica A Juss . J. Basic Microbiol. , 51 , 550 – 556 . Google Scholar Crossref Search ADS PubMed WorldCat Winkelmann , G. ( 1991 ) Structural and stereochemical aspects of iron transport in fungi . Biotechnol. Adv. , 8 , 207 – 231 . Google Scholar Crossref Search ADS WorldCat Winkelmann , G. ( 2002 ) Microbial siderophores-mediated transport . Biochem. Soc. Trans. , 30 , 691 – 695 . Google Scholar Crossref Search ADS PubMed WorldCat Winkelmann , G. ( 2007 ) Ecology of siderophores with specia reference to the fungi . Biometals , 20 , 379 – 392 . Google Scholar Crossref Search ADS PubMed WorldCat Zimmermann , M. B. & Hurrell , R. F. ( 2007 ) Nutritional iron deficiency . Lancet , 370 , 511 – 520 . Google Scholar Crossref Search ADS PubMed WorldCat Appendix. A. Solution of equation (3.8) Here, we consider the model equations $$\begin{equation} \frac{\partial C}{\partial t} =\frac{D}{r} \frac{\partial}{\partial r} \left( r \frac{\partial C}{\partial r}\right), \end{equation}$$(A.1a) $$\begin{equation} \frac{\partial V}{\partial t} =\frac{1}{r} \frac{\partial}{\partial r} \left( r \frac{\partial V}{\partial r}\right), \end{equation}$$(A.1b) for |$1<r<R$| with boundary conditions and initial data $$\begin{equation}\hskip-60pt D\frac{\partial C}{\partial r}(1,t) = - 1, \qquad V(1,t)=0, \end{equation}$$(A.2a) $$\begin{equation} C(R,t)=0, \qquad D \frac{\partial C}{\partial r}(R,t) = -\frac{\partial V}{\partial r}(R,t), \end{equation}$$(A.2b) $$\begin{equation}\hskip-80pt C(r,0)=V(r,0)=0. \end{equation}$$(A.2c) Due to the boundary conditions (A.2b), we first solve equation (A.1a) and then construct the solution for equation (A.1b). A.1 Solution of (A.1a) From the non-homogeneous boundary conditions in equation ( A.2), we suppose that |$C(r,t) = C_{S}(r) + \hat{C}(r,t)$| where |$C_S(r)$| denotes the steady-state distribution and satisfies those same non-homogeneous boundary conditions while |$\hat{C}(r,t)$| satisfies homogeneous boundary conditions (and therefore represents the transition of the initial data (A.2c) towards the final steady state distribution |$C_S(r)$|⁠). The steady-state distribution |$C_S(r)$| satisfies $$\begin{align*} 0 = \frac{\mbox{d}}{\mbox{d} r} \left(r \frac{\mbox{d} C_S}{\mbox{d} r} \right), \qquad \qquad \mbox{for}\ 1<r<R. \end{align*}$$ After integrating twice and applying the non-homogeneous boundary conditions ( A.2), the corresponding steady-state solution is given by $$\begin{equation} C_S(r) = \frac{1}{D} \ln\left(\frac{R}{r}\right). \end{equation}$$(A.3) The function |$\hat{C}(r,t)$| satisfies equation (A.1a) but the corresponding boundary conditions (A.2a) and (A.2b) are expressed as $$\begin{equation} \frac{\partial \hat{C}}{\partial r}(1,t)=0, \qquad \qquad \hat{C}(R,t)=0, \end{equation}$$(A.4) while the corresponding initial data is $$\begin{equation} \hat{C}(r,0) = - C_S(r). \end{equation}$$(A.5) By assuming |$\hat{C}(r,t)=\hat{F}(r)\hat{G}(t)$|⁠, separating variables yields $$\begin{equation} \hat{G}(t) = e^{-\lambda D t} \end{equation}$$(A.6) and the Bessel differential equation |$r \hat{F}^{\prime \prime } + \hat{F}^{\prime } + \lambda r \hat{F}=0$| (where |$^\prime$| denotes differentiation with respect to |$r$|⁠), which has general solution $$\begin{equation} \hat{F} = c_1 J_0\left( \sqrt{\lambda} r\right) + c_2 Y_0\left(\sqrt{\lambda} r\right), \end{equation}$$(A.7) where |$J_0$| and |$Y_0$| are the Bessel functions of first and second kind, respectively. From the boundary condition at |$r=R$|⁠, equation (A.4) allows the constant |$c_2$| to be expressed in terms of |$c_1$| and by introducing |$A=c_2 Y_0(\sqrt{\lambda } R)$| and $$\begin{equation} \phi(r) = \frac{Y_0\left(\sqrt{\lambda}r\right)}{Y_0\left(\sqrt{\lambda}R\right)} - \frac{J_0\left(\sqrt{\lambda} r\right)}{J_0\left(\sqrt{\lambda} R\right)}, \end{equation}$$(A.8) it follows that $$\begin{equation} \hat{F}(r)=A\phi(r). \end{equation}$$(A9) Substituting equations (A.6) and (A.9) into the boundary condition on |$r=1$| yields the characteristic equation $$\begin{equation} \frac{J_1\left(\sqrt{\lambda}\right)}{J_0\left(\sqrt{\lambda}R\right)} - \frac{Y_1\left(\sqrt{\lambda}\right)}{Y_0\left(\sqrt{\lambda} R\right)} =0, \end{equation}$$(A. 10) with roots |$\lambda _n$| denoting the eigenvalues. Thus, |$\hat{C}(r,t)$| can be represented by the series $$\begin{equation} \hat{C}(r,t) = \sum_{n=1}^{\infty} A_n \phi_n(r) e^{-\lambda_n D t}, \qquad \qquad \mbox{for}\ 1<r<R, \end{equation}$$(A.11) where |$A_n$| are constants and the eigenfunctions |$\phi _n(r)$| are obtained from equation (A.8) evaluated with |$\lambda =\lambda _n$|⁠. Notice that |$\phi _n(R)=0$|⁠, |$\phi _n^\prime (1)=0$| and |$r\phi _n^{\prime \prime }+\phi _n^\prime + \lambda _n r \phi _n=0$|⁠. From equation (A.5), there is a generalized Fourier series satisfying $$\begin{align*} \frac{1}{D} \ln \left( \frac{r}{R}\right) = \sum_{n=1}^{\infty} A_n \phi_n(r) \end{align*}$$ and hence the constants |$A_n$| can be determined as $$\begin{equation*} A_n = \frac{\int_{1}^{R} r \ln(r/R) \phi_n(r) \, \mbox{d}r}{ D \int_{1}^{R} r \phi_n(r)^2 \, \mbox{d}r}. \end{equation*}$$ By noting that |$(r \phi _n^\prime )^\prime + \lambda _n r \phi _n=0$|⁠, integrating by parts and recalling that |$\phi _n^\prime (1)=0$| and |$\phi _n(R)=0$|⁠, we obtain $$\begin{equation*} \int_{1}^{R} r \ln(r/R) \phi_n(r) \, \mbox{d}r =-\frac{\phi_n(1)}{\lambda_n}. \end{equation*}$$ Using a similar approach (Bowman, 2003), $$\begin{equation*} \int_{1}^{R} r \phi_n(r)^2 \, \mbox{d}r = \frac{R^2}{2\lambda_n}\left( \phi^\prime_n(R) \right)^2 - \frac{1}{2} \phi_n^2(1). \end{equation*}$$ Consequently, after some simplification, $$\begin{equation} A_n = \frac{2 \phi_n(1)}{D \left[\phi_n^2(1)\lambda_n - R^2 (\phi^\prime_n(R))^2\right]}, \end{equation}$$(A.12) and hence the solution of equation (A.1a) is $$\begin{equation} C(r,t) = \frac{1}{D} \ln \left(\frac{R}{r}\right) + \sum_{n=1}^{\infty} A_n \phi_n(r) e^{-\lambda_n D t}, \qquad \qquad \mbox{for}\ 1<r<R, \end{equation}$$(A13) where |$A_n$|⁠, |$\lambda _n$| and |$\phi _n(r)$| are defined in (A.12), (A.10) and (A.8), respectively. A.2 Solution of (A.1b) Due to the boundary conditions ( A.2), and in particular the flux condition on |$r=R$|⁠, we seek a solution of the form |$V(r,t) = V_S(r) + \bar{V}(r,t) + \hat{V}(r,t)$| where |$V_S(r)$| denotes the steady-state solution, |$\bar{V}(r,t)$| matches the temporal change due to the relationship between the fluxes of |$C(r,t)$| and |$V(r,t)$| at |$r=R$| and |$\hat{V}(r,t)$| accounts for the change from the initial data. The steady-state solution |$V_S(r)$| satisfies the ordinary differential equation $$\begin{equation*} 0 = \frac{\mbox{d}}{ \mbox{d} r} \left( r \frac{\mbox{d} V_S}{\mbox{d} r}\right), \qquad \qquad \mbox{for}\ 1<r<R, \end{equation*}$$ with boundary conditions |$V_S(1)=0$| and |$\frac{\mbox{d} V_S}{\mbox{d} r}(R)=-D \frac{\mbox{d} C_S}{\mbox{d} r}(R)$|⁠. Consequently, we see that $$\begin{equation} V_S(r) = \ln\left( r\right). \end{equation}$$(A.14) The function |$\bar{V}(r,t)$| satisfies the PDE in (A.1b) but with boundary conditions $$\begin{equation} \bar{V}(1,t)=0, \qquad \qquad \frac{\partial \bar{V}}{\partial r}(R,t)=-D \frac{\partial \hat{C}}{\partial r}(R,t), \end{equation}$$(A.15) where |$\hat{C}(r,t)$| is defined in equation (A.11). Due to the form of |$\hat{C}(r,t)$|⁠, suppose $$\begin{align} \bar{V}(r,t) = \sum_{n=1}^{\infty} B_n \psi_n(r) e^{-\lambda_n D t} \end{align}$$(A.16) for suitable constants |$B_n$| and eigenfunction |$\psi _n(r)$|⁠. Since |$\bar{V}(r,t)$| satisfies equation (A.1b), it follows that $$\begin{equation*} r \psi^{\prime\prime} + \psi_n^\prime + \lambda_n D r \psi_n=0 \end{equation*}$$ and hence $$\begin{equation*} \psi_n(r) = c_3 J_0\left( \sqrt{\lambda_n D} r \right) + c_4 Y_0\left(\sqrt{\lambda_n D}r \right), \end{equation*}$$ where |$c_3$| and |$c_4$| are constants. Since |$\psi _n(1)=0$| from (A.15), the constant |$c_3$| can be expressed in terms of |$c_4$| and by substituting into the boundary condition on |$r=R$| and defining $$\begin{equation} \psi_n(r) = \frac{Y_0 \left( \sqrt{\lambda_n D} r \right)}{Y_0 \left( \sqrt{\lambda_n D} \right)} - \frac{J_0 \left(\sqrt{\lambda_n D} r \right)}{J_0 \left( \sqrt{\lambda_n D} \right)}, \end{equation}$$(A.17) it follows that $$\begin{equation} B_n \psi_n^\prime(R) = - D A_n \phi_n^\prime(R) \end{equation}$$(A.18) and provided |$\psi _n^\prime (R) \neq 0$| the constants |$B_n$| can be evaluated. Hence, $$\begin{equation} \bar{V}(r,t) = -D \sum_{n=1}^{\infty} A_n \frac{\phi_n^\prime(R)}{\psi_n^\prime(R)} \psi_n(r) e^{-\lambda_n D t}, \qquad \qquad \mbox{for}\ 1<r<R. \end{equation}$$(A.19) The function |$\hat{V}(r,t)$| satisfies equation (A.1b) with homogeneous boundary conditions and initial data given by $$\begin{equation} \hat{V}(1,t)=0, \qquad \qquad \frac{\partial \hat{V}}{\partial r}(R,t)=0, \qquad \qquad \hat{V}(r,0) = - V_S(r) - \bar{V}(r,0). \end{equation}$$(A.20) By supposing |$\hat{V}(r,t)=\tilde{F}(r)\tilde{G}(t)$|⁠, separating variables and integrating gives $$\begin{equation} \tilde{G}(t) = e^{-\mu t}, \end{equation}$$(A.21) while |$\tilde{F}(r)$| satisfies $$\begin{equation*} r \tilde{F}^{\prime\prime} + \tilde{F}^\prime + \mu r \tilde{F}=0. \end{equation*}$$ As above, the general solution for |$\tilde{F}(r)$| can be expressed in terms of Bessel functions while the boundary condition (A.20) on |$r=1$| gives |$\tilde{F}(r)=E\omega (r)$| where |$E$| is a constant and $$\begin{equation} \omega(r) = \frac{Y_0(\sqrt{\mu}r)}{Y_0(\sqrt{\mu} )} - \frac{J_0(\sqrt{\mu}r)}{J_0(\sqrt{\mu} )}. \end{equation}$$(A.22) The boundary condition (A.20) on |$r=R$| therefore gives the eigenvalues |$\mu _n$| as the roots of $$\begin{equation} \frac{J_1(\sqrt{\mu} R)}{J_0(\sqrt{\mu} )} - \frac{Y_1(\sqrt{\mu} R)}{Y_0(\sqrt{\mu})} =0. \end{equation}$$(A.23) Hence, |$\hat{V}(r,t)$| is given by $$\begin{equation} \hat{V}(r,t) = \sum_{n=1}^{\infty} E_n \omega_n(r) e^{-\mu_n t}, \qquad \qquad \mbox{for}\ 1<r<R, \end{equation}$$(A.24) where |$E_n$| are constants, and |$\omega _n(r)$| is equation (A.22) evaluated at |$\mu =\mu _n$|⁠. Note that for all |$n$|⁠, |$\omega _n(1)=0$|⁠, |$\omega _n^\prime (R)=0$| and |$r\omega _n^{\prime \prime } + \omega _n^\prime + \mu _n r \omega _n=0$|⁠. It now only remains to determine the constants |$E_n$|⁠. The initial data in (A.20) gives $$\begin{equation*} - \ln\left(r\right) + D \sum_{n=1}^{\infty} A_n \frac{\phi^\prime_n(R)}{\psi_n^\prime(R)} \psi_n(r) = \sum_{m=1}^{\infty} E_m \omega_m(r). \end{equation*}$$ By noting that $$\begin{equation*} \int_{1}^{R} r \omega_n(r) \omega_m(r) \, \mbox{d}r =0\ \mbox{ for all}\ \ n\neq m, \end{equation*}$$ it follows that $$\begin{equation} E_m = \frac{\int_{1}^{R} r \left[ - \ln (r) + D \sum_{n=1}^{\infty} A_n \frac{\phi_n^\prime(R)}{\psi_n^\prime(R)}\psi_n(r) \right] \omega_m(r) \, \mbox{d}r}{\int_{1}^{R} r \omega_m^2(r) \, \mbox{d}r}. \end{equation}$$(A.25) By making use of the identities noted above, integration by parts yields $$\begin{equation*}\begin{split} \int_{1}^{R} r \ln (r) \omega_m(r) \, \mbox{d}r &= \frac{\omega_m(R)}{\mu_m},\\ \int_{1}^{R} r \omega_m^2(r) \, \mbox{d}r &= \frac{R^2}{2} \omega_m^2(R) - \frac{1}{2\mu_m} \left(\omega_m^\prime(1)\right)^2,\\ \int_{1}^{R} r \psi_n(r) \omega_m(r) \, \mbox{d}r &=\frac{R \psi_n^\prime(R) \omega_m(R)}{\mu_m - D \lambda_n}, \end{split} \end{equation*}$$ provided |$\mu _m \!\neq\! D \lambda _n$|⁠. By using the above integrals and after some lengthy algebra, equation (A.25) yields $$\begin{equation} E_m = \frac{2 \omega_m(R) \left[ 1 - R \mu_m D \sum_{n=1}^{\infty} \frac{A_n \phi_n^\prime(R)}{\mu_m - D \lambda_n} \right]}{ \left(\omega_m^\prime(1)\right)^2 - \mu_m R^2 \omega_m^2(R)}, \end{equation}$$(A.26) again provided |$\mu _m \neq D \lambda _n$|⁠. Finally from equations (A.14), (A.19) and (A.24), it follows that the solution of equation (A.1b) is given by $$\begin{equation} V(r,t) = \ln\left(r\right) -D \sum_{n=1}^{\infty} A_n \frac{\phi_n^\prime(R)}{\psi_n^\prime(R)} \psi_n(r) e^{-\lambda_n D t} + \sum_{n=1}^{\infty} E_n \omega_n(r) e^{-\mu_n t}, \qquad \qquad \mbox{for}\ 1<r<R, \end{equation}$$(A.27) where |$\phi _n(r)$| is given in equation (A.8), |$\psi _n(r)$| is given in equation (A.17), |$\lambda _n$| are the roots of (A.10), |$\omega _n(r)$| is stated in (A.22), |$\mu _n$| are the roots of (A.23) and the constants |$A_n$| and |$E_n$| are defined by equations (A.12) and (A.26), respectively. B. Derivation of approximations (3.15) The characteristic equations (3.11) are of the form $$\begin{equation} J_1(x p) Y_0(x q) - J_0(x q) Y_1(x p)=0, \end{equation}$$(B.1) where |$x$| denotes the square root of the eigenvalue and |$p\neq q$| take either the values |$1$| or |$R$| (depending on which eigenvalue |$\lambda$| or |$\mu$| is being considered). Using Hankel’s asymptotic expansions, Harrison (2009) obtained approximations, valid for large values of |$z$|⁠, for the Bessel functions $$\begin{eqnarray*} J_n(z) & = & \sqrt{\frac{2}{\pi z}} \beta_n(z) \cos\left(z - \frac{\pi}{4} - \alpha_n(z) \right), \\ Y_n(z) & = & \sqrt{\frac{2}{\pi z}} \beta_n(z) \sin\left(z - \frac{\pi}{4} - \alpha_n(z) \right), \end{eqnarray*}$$ for suitable series |$\alpha _n(z)$| and |$\beta _n(z)$| each in terms of powers of |$1/z$|⁠. To determine the roots of the characteristic equation (B.1), we note from Harrison (2009) that for |$|z|\gg 1$| $$\begin{eqnarray*} \alpha_0(z) & = & \frac{1}{8z} - \frac{25}{384z^3} + \frac{1073}{5120 z^5} - \frac{375733}{229376 z^7} + O\left(\frac{1}{z^9}\right), \\ \alpha_1(z) & = & -\frac{3}{8z} + \frac{21}{128 z^3} - \frac{1899}{5120 z^5} + \frac{543483}{229376 z^7} + O\left(\frac{1}{z^9}\right). \end{eqnarray*}$$ By using these expressions in the generalized characteristic equation (B.1), it follows that $$\begin{equation} \sin \left( (q-p)x + \frac{\pi}{2} + \alpha_1(px) - \alpha_0(qx) \right) = 0. \end{equation}$$(B.2) Consequently, equation (B.2) gives rise to $$\begin{equation} x \left[ 1+ \frac{\alpha_1(px) - \alpha_0(qx)}{(q-p)x} \right] = P_n(p,q), \end{equation}$$(B.3) where |$n$| is an integer and $$\begin{equation*} P_n(p,q) = \frac{\pi (n-\frac{1}{2})}{q-p}. \end{equation*}$$ Since we seek positive roots |$x$| and note that $$\begin{equation*} \frac{\alpha_1(px)-\alpha_0(qx)}{(q-p)x} \to 0 \qquad \mbox{as}\ x\to \infty, \end{equation*}$$ it follows that |$P_n(p,q)>0$| and so if |$q>p$| then |$n$| has to be a positive integer. (The case |$q<p$| is treated below.) Equation (B.3) can therefore be written as the summation of a series of even powers of |$1/x$| $$\begin{equation*} x\left[ 1+ \frac{a_2(p,q)}{x^2} + \frac{a_4(p,q)}{x^4} + O\left(\frac{1}{x^6}\right) \right] = P_n(p,q), \end{equation*}$$ where the coefficients |$a_n(p,q)$| are easily determined from the above series for |$\alpha _0(qx)$| and |$\alpha _1(px)$|⁠. By constructing the reciprocal of the series on the left, it follows that $$\begin{equation*} \begin{split} &\frac{1}{P_n(p,q)} = \frac{1}{x} + \frac{3q+p}{8pq(q-p) x^3} + \frac{25p^4 - 19p^3q+36p^2q^2+117pq^3-63q^4}{384 p^3 q^3 (q-p)^2 x^5} \\ & \; + \frac{3219p^7\!-6188p^6q+\!3749p^5q^2\!-\!480p^4q^3\!+\!1440p^3q^4\!+\!7767p^2q^5\!-\!13284pq^6\!+5697q^7}{15360 p^5 q^5 (q-p)^3 x^7} \!+\! O\left(\frac{1}{x^9}\right). \end{split} \end{equation*}$$ By using series inversion, a power series for |$\frac{1}{x}$| in terms of odd powers of |$\frac{1}{P_n(p,q)}$| is obtained and is given by $$\begin{equation*} \begin{split} \frac{1}{x} =&\ \frac{1}{P_n(p,q)} - \frac{p+3q}{8pq(q-p) P_n(p,q)^3} - \frac{25p^4 - 37p^3q -72p^2q^2-45pq^3-63q^4}{384 p^3q^3(q-p)^2 P_n(p,q)^5} \\ & \; - \frac{1073p^7 - 2396p^6 q+623p^5q^2 +1200p^4q^3+720p^3q^4+1989p^2q^5-1908pq^6+1899q^7}{5120p^5q^5(q-p)^3 P_n(p,q)^7}\\ &+O\left(\frac{1}{P_n(p,q)^9}\right). \end{split} \end{equation*}$$ Multiplying through by |$P_n(p,q)$| produces a power series for |$\frac{P_n(p,q)}{x}$| in terms of even powers of |$\frac{1}{P_n(p,q)}$|⁠. Taking reciprocals and then finally multiplying through by |$P_n(p,q)$| yields $$\begin{equation*} x = P_n(p,q) + \frac{Q_1(p,q)}{P_n(p,q)} + \frac{Q_3(p,q)}{P_n(p,q)^3} + \frac{Q_5(p,q)}{P_n(p,q)^5} + O\left( \frac{1}{P_n(p,q)^7}\right), \end{equation*}$$ where |$P_n(p,q)$| is as defined above and $$\begin{equation*} \begin{split} Q_1(p,q) &= \frac{p+3q}{8 p q (q-p)}, \\ Q_3(p,q) & = \frac{25 p^4-31 p^3 q-36 p^2 q^2+9 p q^3-63 q^4}{384 (q-p)^2 p^3 q^3}, \\ Q_5(p,q) & = \frac{3219p^7 - 6938p^6q + 2279p^5 q^2 + 2040p^4q^3+360p^3q^4+4797p^2q^5-7614pq^6+5697q^7}{15360 p^5 q^5 (q-p)^3}. \end{split} \end{equation*}$$ If on the other hand, |$q<p$| then from (B.2) $$\begin{equation*} x \left[ 1+ \frac{\alpha_1(px) - \alpha_0(qx)}{(q-p)x} \right] = \bar{P}_n(p,q), \end{equation*}$$ where $$\begin{equation*} \bar{P}_n(p,q) = \frac{\pi (\frac{1}{2}-n)}{q-p} \end{equation*}$$ and |$n$| is now a strictly positive integer. Notice that |$\bar{P}_n(p,q) = P_{n}(q,p)$| and the remaining terms in the expansion are obtained in the same way as for the case |$q>p$|⁠. C. Derivation of (3.16) and (3.17) Here, approximations for eigenvalues |$\lambda _n$| and |$\mu _n$| in equation (3.11) are derived for the case of large |$R$|⁠. Attention is focussed on the smallest eigenvalues since they exert the greatest influence on the solution (3.10). The following Bessel function expansions, valid as |$x \to 0$|⁠, will be used $$\begin{equation} \begin{split} J_0(x) & = 1 -\frac{x^2}{4} + O(x^4), \\ Y_0(x) & = \frac{2}{\pi} \ln\left( \frac{x e^{\gamma}}{2} \right) -\frac{x^2}{2\pi} \left[ \ln\left( \frac{x e^{\gamma}}{2} \right) -1 \right] +O(x^4\ln(x)), \\ J_1(x) &= \frac{x}{2} -\frac{x^3}{16} +\frac{x^5}{384} + O(x^7), \\ Y_1(x) &= -\frac{2}{\pi x} + \frac{x}{\pi} \left[ \ln\left( \frac{x e^{\gamma}}{2} \right) -\frac{1}{2} \right] \!- \frac{x^3}{ 8\pi} \left[ \ln\left( \frac{x e^{\gamma}}{2} \right) \!-\frac{5}{4} \right] + \frac{x^5}{192\pi} \left[ \ln\left( \frac{x e^{\gamma}}{2} \right) -\frac{5}{3} \right] \!+ O(x^7\ln(x)), \end{split} \end{equation}$$(C.1) where |$\gamma$| denotes Euler’s constant. C.1 Derivation of (3.16): approximation for small |$\lambda _n$| We recall that |$\lambda _n$| satisfies $$\begin{equation} J_1 \left(\sqrt{\lambda_n}\right) Y_0 \left(R \sqrt{\lambda_n}\right) = Y_1\left(\sqrt{\lambda_n}\right) J_0\left(R \sqrt{\lambda_n}\right). \end{equation}$$(C.2) As |$R \rightarrow \infty$|⁠, by numerically solving equation (C.2), we find that |$\lambda _n \rightarrow 0$|⁠. First, we expand the two functions of only |$\sqrt{\lambda _n}$| using the expansions for |$J_1(x)$| and |$Y_1(x)$| in (C.1). Then multiplying by |$2\pi \sqrt{\lambda _n}$| yields $$\begin{equation} \left( \pi \lambda_n + O({\lambda_n}^2) \right) Y_0 \left(R \sqrt{\lambda_n}\right) = \left( -4 + O(\lambda_n\ln(\lambda_n)) \right) J_0 \left(R \sqrt{\lambda_n}\right). \end{equation}$$(C.3) By numerically solving equation (C.2), we find that |$R \sqrt{\lambda _n}$| tends to a constant as |$R \rightarrow \infty$|⁠, so we seek an expansion in the form $$\begin{equation} R \sqrt{\lambda_n} = \zeta_n + \epsilon. \end{equation}$$(C.4) We use the Taylor series expansions as |$\epsilon \rightarrow 0$| and substituting into (C.3) yields $$\begin{equation*} \begin{split} \left( \pi \frac{\zeta_n^2} { R^2} + O(\epsilon R^{-2}, R^{-4}) \right) & \left( Y_0(\zeta_n) +\epsilon Y_0^{\prime}(\zeta_n) + O(\epsilon^2) \right) \\ & = \left( -4 + O(R^{-2}\ln(R)) \right) \left( J_0(\zeta_n) +\epsilon J_0^\prime(\zeta_n) + O(\epsilon^2) \right). \end{split} \end{equation*}$$ Notice that in the equation above, as |$R \rightarrow \infty$| the left-hand side tends to zero but the right-hand side tends to |$-4J_0(\zeta _n)$|⁠. Thus, we require $$\begin{equation} J_0(\zeta_n)=0. \end{equation}$$(C.5) Hence, the |$\zeta _n$|’s in (C.4) are the |$n^{\tiny{th}}$| roots of |$J_0$|⁠. Using this and keeping the leading order terms yields $$\begin{align*} \frac{\pi \zeta_n^2} { R^2} Y_0(\zeta_n) + O(\epsilon R^{-2}, R^{-4}) = -4 \epsilon J_0^{\prime}(\zeta_n) + O( \epsilon R^{-2} \ln(R), \epsilon^2). \end{align*}$$ Thus, to leading order, $$\begin{equation*} \epsilon = -\frac{\pi \zeta_n^2 Y_0(\zeta_n)} { 4R^2 J_0^\prime(\zeta_n)} +O(R^{-4}\ln(R)). \end{equation*}$$ Hence, from (C.4) it follows that $$\begin{equation} \sqrt{\lambda_n} = \frac{\zeta_n} { R} -\frac{\pi \zeta_n^2 Y_0(\zeta_n) }{ 4R^3 J_0^\prime(\zeta_n)} + O(R^{-5}\ln(R)), \end{equation}$$(C.6) which is valid for |$\sqrt{\lambda _n} \ll 1$|⁠, i.e. |$\zeta _n \ll R$|⁠. C.2 Approximation for small |$\mu _n$|⁠, |$n \ge 2$| Recall that the eigenvalue |$\mu _n$| satisfies $$\begin{equation} J_1 \left(R \sqrt{\mu_n}\right) Y_0\left(\sqrt{\mu_n}\right) = Y_1\left(R \sqrt{\mu_n}\right) J_0\left(\sqrt{\mu_n}\right). \end{equation}$$(C.7) As |$R \rightarrow \infty$|⁠, by numerically solving equation (C.7), we find that |$\mu _n \rightarrow 0$|⁠. First, we expand the two functions of only |$\sqrt{\mu _n}$| using the series for |$J_0(x)$| and |$Y_0(x)$| in equation (C.1) so that equation (C.7) becomes $$\begin{equation} J_1(R \sqrt{\mu_n}) \left( \frac{2} { \pi} \ln\left( \frac{\sqrt{\mu_n} e^{\gamma}}{ 2} \right) +O(\mu_n\ln(\mu_n)) \right) = Y_1(R \sqrt{\mu_n}) (1 + O(\mu_n)). \end{equation}$$(C.8) By numerically solving equation (C.7), we find that |$R \sqrt{\mu _n}$| tends to a constant as |$R \rightarrow \infty$|⁠, so we seek an expansion in the form $$\begin{equation} R \sqrt{\mu_n} = \theta_{n} + \delta. \end{equation}$$(C.9) A Taylor series expansion as |$\delta \to 0$| is constructed from equation (C.8) resulting in $$\begin{equation*} \begin{split} \left( J_1( \theta_{n}) + \delta J_1^\prime(\theta_{n}) + O(\delta^2) \right) & \left( \frac{2} { \pi} \ln\left( \frac{(\theta_{n} + \delta) e^{\gamma}}{ 2R} \right) +O( R^{-2} \ln(R) ) \right) \\ &= \left( Y_1(\theta_{n}) + \delta Y_1^\prime(\theta_{n}) + O(\delta^2) \right) \left( 1 + O\left( R^{-2} \right) \right). \end{split} \end{equation*}$$ We notice that in the equation above, as |$R \rightarrow \infty$| the right-hand side remains finite but the left-hand side tends to infinity like |$-2\ln (R)J_1(\theta _{n})/\pi$|⁠. Thus, we require $$\begin{equation} J_1(\theta_{n})=0. \end{equation}$$(C.10) Hence, the |$\theta _n$|’s in (C.9) are the |$n^{\tiny{th}}$| roots of |$J_1$| and note that |$\theta _1=0$| is the first solution. Before collecting leading order terms, notice that the approach fails around |$\theta _1$| since |$Y_1(0)$| is undefined and hence an alternative approach is required for the calculation of |$\mu _1$| (see subsection C.3). Provided |$n \geq 2$|⁠, keeping the leading order terms yields $$\begin{equation*} -\delta J_1^\prime(\theta_{n}) \frac{2} { \pi} \ln(R) = Y_1(\theta_{n}) + O\left( \delta, R^{-2}, R^{-2} \ln(R) \delta \right). \end{equation*}$$ Thus, to leading order and provided |$1 \ll \ln (R)$|⁠, i.e. |$e \ll R$|⁠, $$\begin{align*} \delta = -\frac{\pi Y_1(\theta_{n})} { 2 J_1^\prime(\theta_{n}) \ln(R)} + O\left( \frac{1 }{ \ln(R)^2}, \frac{R^{-2}}{ \ln(R)} \right). \end{align*}$$ Hence from (C.9), we have $$\begin{equation} \sqrt{\mu_n} = \frac{\theta_{n}} { R} -\frac{\pi Y_1(\theta_{n}) }{ 2R \ln(R) J_1^\prime(\theta_{n})} + O\left( \frac{R^{-1}}{ \ln(R)^2}, \frac{R^{-3}}{ \ln(R)} \right), \end{equation}$$(C.11) which is valid for |$\sqrt{\mu _n} \ll 1$|⁠, i.e. |$\theta _{n} \ll R$| (and the condition |$e\ll R$| is ensured since |$e<\theta _2$|⁠). C.3 Derivation of (3.17): approximation for small |$\mu _1$| The above approach failed to calculate |$\mu _1$| because |$Y_1(0)$| is not defined and hence an alternative approach, utilizing a different expansion, is described here. Recall |$\mu _1$| satisfies $$\begin{equation} J_1(R \sqrt{\mu_1}) Y_0(\sqrt{\mu_1}) = Y_1(R \sqrt{\mu_1}) J_0(\sqrt{\mu_1}). \end{equation}$$(C.12) As |$R \rightarrow \infty$|⁠, by numerically solving equation (C.12), we find that |$R \sqrt{\mu _1} \rightarrow 0$|⁠. By substituting all the expansions in equation (C.1) into equation (C.12), it follows that $$\begin{equation*} \begin{split} &\left( {R \sqrt{\mu_1} \over 2} -{R^3 \mu_1^{3 \over 2} \over 16} +{R^5 \mu_1^{5 \over 2} \over 384} + O(R^7 \mu_1^{7 \over 2}) \right) \left( {2 \over \pi} \ln\left( {\sqrt{\mu_1} e^{\gamma} \over 2} \right) + O( \mu_1 \ln(\mu_1) ) \right) \\ & \qquad = \left( 1 + O(\mu_1) \right) \times \left( \!\! -{2 \over \pi R \sqrt{\mu_1}} \!+\! {R \sqrt{\mu_1} \over \pi} \! \left[ \ln\left( \! {R \sqrt{\mu_1} e^{\gamma} \over 2} \right) \!-\! {1 \over 2} \right] \right. \\ & \qquad \qquad \left. \!- {R^3 \mu_1^{3 \over 2} \over 8\pi} \! \left[ \ln\left( \! {R \sqrt{\mu_1} e^{\gamma} \over 2} \! \right) \!-\! {5 \over 4} \right] \!+\! {R^5 \mu_1^{5 \over 2} \over 192\pi} \! \left[ \ln\left( \! {R \sqrt{\mu_1} e^{\gamma} \over 2} \! \right) \!-\! {5 \over 3} \right] \!+\! O(R^7 \mu_1^{7\over 2} \ln(R \sqrt{\mu_1})) \!\! \right). \end{split} \end{equation*}$$ Multiplying by |$2\pi R \sqrt{\mu _1}$| and cancelling out terms reduces this expression to $$\begin{equation} 0= -4 + R^2 \mu_1 \left[ 2\ln(R) - 1 \right] - {R^4 \mu_1^2 \over 4} \left[ \ln(R) - {5 \over 4} \right] + {R^6 \mu_1^3 \over 96} \left[ \ln(R) - {5 \over 3} \right] + O(\mu_1, R^8 \mu_1^4 \ln(R \sqrt{\mu_1})). \end{equation}$$(C.13) Next, motivated by the presence of |$\ln (R)$| and the powers of |$\mu _1$| in the above, we suppose that |$\mu _1$| can be expanded in the form $$\begin{equation} \mu_1 = {1 \over R^2 \ln(R)} \left[ a + {b \over \ln(R)} + {c \over \ln(R)^2} + O\left( {1 \over \ln(R)^3} \right) \right], \end{equation}$$(C.14) where |$a,b$| and |$c$| are constants to be determined. By substituting equation (C.14) into equation (C.13) and retaining leading order terms yields $$\begin{align*} O\left( {1 \over R^2 \ln(R)}, {\ln(\ln(R)) \over \ln(R)^3} \right) =&\ -4 \!+\! 2a \!+\! {2b \over \ln(R)} \!+\! {2c \over \ln(R)^2} \!-\! {a \over \ln(R)} \!-\! {b \over \ln(R)^2} \\ & -\! {a^2 \over 4 \ln(R)} \!-\! {ab \over 2 \ln(R)^2} \!+\! {5 a^2\over 16 \ln(R)^2} \!+\! {a^3 \over 96 \ln(R)^2}. \end{align*}$$ Finally, by equating the coefficients of the powers of |$\ln (R)$|⁠, values for |$a,b$| and |$c$| can be determined and hence $$\begin{equation} \mu_1 = {1 \over R^2 \ln(R)} \left[ 2 + {3 \over 2\ln(R)} + {5 \over 6\ln(R)^2} + O\left( {1 \over \ln(R)^3} \right) \right], \end{equation}$$(C.15) which is valid for |$e\ll R$|⁠. © The Author(s) 2020. Published by Oxford University Press on behalf of the Institute of Mathematics and its Applications. All rights reserved. This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)
TI - Mathematical modelling of fungi-initiated siderophore–iron interactions
JF - Mathematical Medicine And Biology: A Journal Of The Ima
DO - 10.1093/imammb/dqaa008
DA - 2020-12-15
UR - https://www.deepdyve.com/lp/oxford-university-press/mathematical-modelling-of-fungi-initiated-siderophore-iron-05jbfd2sje
SP - 515
EP - 550
VL - 37
IS - 4
DP - DeepDyve
ER -