Microbial mutualism at a distance

Complex microbial communities play essential roles in the proper functioning of the environment, in maintaining the health of plants and animals, and in many industrial processes. Within these communities, microbial interactions are often predicated on metabolite exchange, as auxotrophs depend on nutrients made by other microbes. Here, we investigate a mathematical model of growth and interactions between mutualistic microbial species separated in space but coupled by a channel through which nutrients are exchanged diffusively. The model is used to study mutualistic algal-bacterial interactions, focusing on a synthetic model system. Solutions to the model reveal rich dynamics, and allow prediction of the conditions for the successful establishment of remote mutualisms. We connect our findings to understanding complex behaviour in synthetic and naturally occurring microbial communities.


Introduction
Microorganisms display a broad spectrum of interactions critical to the function and dynamics of microbial communities (Widder et al. 2016). A wealth of experimental data on such dynamics is now available from 'omics' approaches (Cooper & Smith 2015, Widder et al. 2016), but these must be complemented by lab-based studies of synthetic consortia and mathematical modelling to reach a mechanistic understanding of microbial dynamics (Widder et al. 2016, Abreu & Taga 2016. The focus of the present work is on mutualistic interactions within microbial communities, as found in several recent examples across the kingdoms of life. These include strains of enteric bacteria (Harcombe et al. 2014) and yeast (Allen et al. 2013) that have been engineered to be mutualistic, and synthetic consortia combining wild type microbial species, such as bacterial tricultures (Kim et al. 2008), mixed cultures of algae and fungi (Hom & Murray 2014), and algae and bacteria (Croft et al. 2005, Kazamia et al. 2012, Wang et al. 2014, Segev et al. 2016.
Historically, mathematical studies of mutualism have been scarce, with models concentrating on predator-prey (Lotka-Volterra) or competitive interactions (Murray 1989). As in Lotka-Volterra systems, models of mutualism often take the form of coupled ordinary differential equations describing the population dynamics of two or more species (Murray 1989, Holland & Deangelis 2010, within network that can be engineered in the laboratory. The nodes of this physically structured network represent reservoirs of different volumes filled with different growing microbial species diffusively exchanging metabolites via porous channels, as described in the model formulated in this work. Diffusive exchanges are parameterised by sets of geometric parameters, as such as the lengths, λ ij , of the channels connecting nodes. which mutualistic interactions appear implicitly through a beneficial effect on partner growth parameters [e.g. by making carrying capacity a function of the partner concentration (Murray 1989, Yukalov et al. 2012, Grant et al. 2014]. These implicit models have significant limitations, particularly when the mutualism arises from the exchange of limiting nutrients, for if the nutrient dynamics are not modelled explicitly, the effect of a partner is difficult to parameterise (Momeni et al. 2017). Models of cooperative synthesis in bacteria have considered the mutualistic interactions between different species through explicit nutrient exchanges in well-mixed cultures (Lee et al. 1976, Wang et al. 2007). To investigate spatially structured mutualistic systems (i.e. those that are not well-mixed), spatial transport of nutrients needs to be included explicitly.
Recent studies have considered spatial aspects of mutualistic interactions. Simulations using flux balance analysis (FBA) successfully predict the spatial growth on agar of colonies of synthetically mutualistic enteric bacteria (Harcombe et al. 2014). The FBA approach requires explicit knowledge of every known metabolic biochemical pathway in each mutualistic species, restricting its applicability to synthetic mutualisms between metabolically well-characterised organisms. Spatial effects on cheating (Momeni et al. 2013) and genetic drift (Müller et al. 2014) observed in yeast colonies growing on agarose pads have also been modelled explicitly. In these models, coupled cells and nutrients diffusing in two dimensions are simulated to predict how nutrient-mediated interactions control spatial heterogeneity and survival of the populations. However, the homogenous environment of nutrient agarose does not possess the intrinsic geometric or topological structure of real microbial environments, such as the matrix of soil or a biofilm. Mutualistic microbial dynamics have not thus far been studied in such structured environments.
Here, we study a model of mutualistic microbial species in a simple geometry representing a minimal unit for a structured environment: populations growing in spatially separated reservoirs, metabolically linked by a channel. As soil can be approximated at the microbial scale as a physical network of growth chambers linked by channels (Pérez-Reche et al. 2012), this model provides the basis for the description of complex natural networks or artificial regular networks (see Figure 1) and the microbial interactions therein. The model is generally applicable to auxotrophs cross-feeding remotely. We apply it to make predictions for the dynamics of mutualistic populations of algae and bacteria exchanging vitamin B 12 and a carbon source (Kazamia et al. 2012). Our predictions provide new insights into the behaviour of microbial communities residing in structured geometries, both within synthetic consortia in the laboratory and environmental microbial communities. Figure 2: Diffusive cross-feeding at a distance. Auxotrophic microbial populations A and B (concentrations a andb) reside in well-mixed reservoirs of equal volume Γ separated by a channel of length L and crosssection Σ. Microbe A produces a carbon source C, of homogeneous concentrationc a , in its reservoir. This diffuses through the channel, forming a profilec(x,t), a function of position along the channelx and timet. On reaching the reservoir where microbe B resides the concentration is homogenised toc b . Symmetrically, the vitamin V produced by microbe B in its reservoir at concentrationv b , diffuses to reservoir A creating a profilev(x,t), homogenised tov a in the reservoir. Here, this general model is applied to an algal-bacterial partnership.

Model
The model describes two populations of mutualistic microbial species, A and B, interacting at a distance. The mutualistic interactions are predicated on auxotrophy: A requires metabolite V (for "vitamin"), excreted by B; conversely B requires metabolite C (for "carbon"), excreted by A. In formulating the problem we shall first use variables with an overbar to denote dimensional quantities (concentrations, time, space), reserving symbols without typographical modification for appropriately rescaled variables. Populations of A and B, with densitiesā(t) andb(t) respectively, reside in two well-mixed reservoirs, of equal volume Γ. These are spatially separated, but connected by a cylindrical channel (length L, cross-sectional area Σ), as in Figure 2. The channel is impervious to cells, but porous to metabolite exchange by diffusion. Population A produces metabolite C with local concentrationc a (t), which diffuses out of the reservoir and into the channel atx = 0 (withx denoting the position along the channel axis), where it develops a spatial profilec(x,t) and eventually reaches the other reservoir atx = L, where its concentration isc b (t). Symmetrically, metabolite V produced by B with concentrationv b (t), diffuses out atx = L givingv(x,t), feeding the other reservoir atx = 0, generating a concentrationv a (t).
We first consider dynamics within the channel connecting the reservoirs, within which metabolites obey one-dimensional diffusion equations, with D s the diffusion coefficients for metabolite S = C or V. The boundary conditions to (2.1) obtained from continuity at the channel-reservoir interface are:c a (t) =c(0,t),c b (t) =c(L,t), v a (t) =v(0,t),v b (t) =v(L,t). Clearly, one characteristic time scale of the problem is set by diffusive equilibration along the length of the channel, where we anticipate that the diffusion constants of both metabolite species are similar. From Fick's law, the flux J s (molecules area −1 time −1 ) of metabolite species S (C or V) entering, say, the left reservoir from the channel is The rate such molecules enter the reservoir is J 0 s Σ, and with instantaneous homogenenisation there, the rate of change of the reservoir concentrations a is J 0 s Σ/Γ. The characteristic length will play an important role in the model. If ∆s is a typical difference in concentration of S between the two reservoirs, then the typical gradient within the channel is ∆s/L, giving rise, by the arguments above, to an associated rate of change of reservoir concentration scaling as ds/dt ∼ (Σ/Γ)D s ∆s/L ∼ D s ∆s/ L, from which we can identify a characteristic equilibration time We define the ratio of equilibration and diffusive time scales to be The regime ζ 1 is that of fast establishment of the linear concentration profile in the tube relative to changes of concentrations in the reservoirs, while for ζ ≥ 1 the transients within the channel are on comparable time scales to that for changes in the reservoirs. Semi-analytical solutions to the problem of chemical diffusion between two connected reservoirs further demonstrate the existence of these two regimes and the role of the previously identified timescales (Materials and Methods).
We now turn to the population dynamics within the reservoirs, in which we explicitly assume that algae reside in reservoir A and bacteria in B, and that vitamin B 12 and carbon are exchanged. The dynamics obey the ordinary differential equations is the flux of metabolite S = C or V entering the right reservoir. In equations (2.7a) we model cell growth as logistic, with maximum growth rate µ i and carrying capacity K i for species i = A or B. Growth rates are limited by the abundance of the required metabolites. This is modelled using Monod factors Monod (1949), e.g., for C, µ bc /(K c +c), where K c is the halfsaturation constant (and symmetrically for V). Linear death terms, with mortality rates δ i for i = A or B, ensure exponential negative growth in the absence of the limiting metabolites. Equations (2.7b) describe the dynamics of metabolite C. This is produced by species A in proportion to its concentration with a rate p c , and diffuses out at 0. In the other reservoir, C is taken up by B. The uptake is assumed proportional to the cell growth rate, the proportionality constant is where Y b is the yield coefficient (how much metabolite C results in a given concentration of species B). Equations (2.7c) describe the V dynamics, which are completely symmetric to the C dynamics. Although inspired by bacterial-algal symbiosis, it is clear that the structure of these dynamics is quite broadly applicable to mutualistic systems in general.

Identifying the key model parameters
In order to access the general dynamics of remotely cross-feeding monocultures, we nondimensionalise equations (2.7). Because our focus is on the impact of geometry on the biological processes, we choose a scheme accordingly. First, normalize the bacterial and algal concentrations by their respective carrying capacities, the organic carbon and vitamin concentrations by their respective half-saturation concentrations, rescale time by the bacterial growth rate, and rescale space by the length scale b = D c /µ b of organic carbon diffusion on the time scale of bacterial growth, defining (2.8) The ratios of algal and bacterial growth rates and of their diffusion constants, are two additional parameters. With now three characteristic lengths in the problem (L, , b ) one can form two independent dimensionless ratios. These can be taken to be (2.10) so that the parameter ζ, defined previously in Eq. 2.6, is ζ = λ/η. There are three pairs of parameters remaining which capture the relative strength of cellular death, uptake and production in bacteria and algae respectively. They are: the ratios of death rate to maximum growth rate of bacteria and algae, which define mortality parameters which must be less than 1 for any population increase to occur; and finally, for both carbon and vitamin, the ratios of the typical uptake rate to the typical rate of change define the uptake parameters for both carbon and vitamin, the ratios of the typical production rate to the typical rate of change define the production strengths (2.13) With these rescalings, the dimensionless evolution equations are for c and v on the interval x ∈ [0, λ], ensuring continuity of fluxes and concentrations at the ends of the tube. Numerical solutions to these dynamics were obtained using an explicit improved Euler scheme for time stepping and a central difference scheme for spatial derivatives. Before discussing parameterisation of this model for a particular mutualistic association, we note, as described below, that this dynamical system exhibits a set of fixed points for which the cell concentrations vanish (a = b = 0) at any combination of residual metabolite concentrations. More interestingly, a single positive fixed point exists given by equations (2.22) below under the conditions of equations (2.23) (see Materials and Methods). This positive fixed point corresponds to an established co-culture at a distance.

Parameterisation for specific microbial associations
The results presented below were obtained from numerical studies of the mathematical model with parameter values corresponding to the mutualistic association between Lobomonas rostrata, a B 12requiring green alga, and Mesorhizobium loti, a B 12 -producing soil bacterium (Kazamia et al. 2012). The following procedure was used to obtain parameter values. First, physiologically relevant ranges for each parameter were collected by direct measurement or from the published literature. Then, specific parameters -both nondimensional parameters of the reduced model and dimensional parameters to convert experimental data to nondimensional units-were obtained by minimizing the squared distance between simulated time evolution, obtained through a custom finite difference solver in Python, and experimental results on mixed cultures, while searching within domains of parameter values which contain the physically relevant ones, and validating the fixed-point conditions in equation (2.23). The basin-hopping minimisation procedure gives local optima which capture well the observed dynamics of mixed co-cultures of L. rostrata and M. loti (see Materials and Methods). The range of physiologically relevant parameters used to constrain the search of parameters for the association of M. loti and L. rostrata are presented in table 1, while the fitted parameters, both dimensional and nondimensional, are given in tables 2 and 3.

Feeding on a distant passive source
Before considering the fully coupled system dynamics, we consider the case of a single auxotrophic species B, concentration b, residing in a reservoir initially free of a growth-limiting metabolite coupled by the channel (also initially nutrient-free) to a strong source of the metabolite with initial concentration c 0 a . This source consists of a reservoir filled with limiting metabolite. The long time steady-state for the model is always extinction of B once it has exhausted the remote resource. However, separation of the microbial population from the source modifies the transient population dynamics. Recalling the nondimensional channel length λ = L/ b and equilibration length η = / b , we can define the nondimensional timescales t diff = λ 2 and t eq = λη as the ratios between the typical times of diffusion and of equilibration between reservoirs, and the biological growth timescale τ b = 1/µ b . These ratios gauge the relative rates of diffusion/equilibration and growth. We require t eq and t diff ∼ 1 for diffusion to transport metabolites to species B, stimulating its growth.
We have solved the remotely-fed single microbe limit of the model (see Materials and Methods) to predict the dynamics of the rhizobial bacterium Mesorhizobium loti fed from a remote glycerol carbon source. Figure 3 shows the transient growth dynamics in the regime for which both geometric parameters λ and η impact the dynamics. We first consider the effect of diffusive reservoir equilibration, quantified by η for a fixed channel length λ. For large η, t eq is large: diffusive equilibration in the reservoir is much slower than growth. Thus, the instantaneous flux from the carbon  Link et al. (2008) † this work: see Materials and Methods. a considering glycerol or small sugars such as glucose and sucrose. b obtained considering E. coli and species of rhizobia growing on different sugars. The range of values is quite wide due to the ability of bacteria to tune their affinity constant depending on the environmental conditions (Ferenci 1999). c obtained considering L. rostrata and other B 12 -dependent species. d obtained considering two species belonging to the same family (Chlamydomonadaceae) as L. rostrata, and arabinose molar mass. e obtained considering two B 12 -producing bacterial species, Azobacter vinelandii and Halomonas sp.

Fitted dimensional parameter Symbol Value
Algal carrying capacity K a 3.0 × 10 6 cells mL −1 Bacterial carrying capacity source reservoir to the bacterial reservoir is below what the bacteria need to grow to carrying capacity. As a result, increasing η decreases the value of the peak bacterial concentration (preceding the inevitable decay), as well as delaying the onset of growth ( Figure 3A). Next we fix η and vary λ. Since the diffusive timescale scales like λ 2 , increasing λ progressively delays the onset of bacterial growth ( Figure 3B inset). Large λ values also correspond to weaker carbon source gradients across the tube, and thus a 'slow-release' nutrient flux. Consequently, a less concentrated population can be sustained for longer by the remote source ( Figure 3B). The passive source case we have just considered demonstrates the critical role played by the geometric parameters λ and η in setting the timescale of transients, but also the peak microbial numbers achievable on a finite resource. a obtained considering carbon with diffusivity D c = 5 × 10 −6 cm 2 s −1 as metabolite C and B 12 vitamin with diffusivity D v = 2 × 10 −6 cm 2 s −1 as metabolite V

Remotely cross-feeding populations
Next, we consider auxotrophic populations in separate reservoirs, exchanging limiting metabolites through a connecting channel. As mentioned earlier, we apply the model to an algal-bacterial system, obtaining our parameters from experiments where the phototrophic alga L. rostrata, auxotrophic for vitamin B 12 , is grown in co-culture with the heterotrophic bacterium M. loti. The algal and bacterial populations in their reservoirs have initial concentrations, a 0 and b 0 , respectively. Neither carbon source nor vitamin (the limiting metabolites) are initially present in the reservoirs and channel. The coexistence diagrams in Figure 4 show what values in the initial concentration parameter space give rise to long-term mutualistic coexistence or a population crash due to metabolite deprivation. These are the possible fixed points of our model, which we shall also refer to as model equilibria. Figure  4A displays the boundary between these two regions for different values of the channel length λ for a fixed value of the equilibration length η. In Figure 4B crash-coexistence boundaries are instead shown for different equilibration lengths η at fixed λ. Also shown on both diagrams is the membrane limit (grey dashed line). In this limit the distance between reservoirs vanishes (λ → 0) and they are simply separated by a membrane impervious to cells, as has been demonstrated experimentally in co-culturing/metabolomic experiments (Paul et al. 2013). We see that increasing the channel length has the effect of pushing the crash-coexistence boundary toward higher initial microbial concentrations. Coexistence is achieved in the membrane limit for initial concentrations lower than those for finite λ. The boundary between crash and coexistence regions shifts quantitatively with λ, but does not change significantly qualitatively. Its shape is revealing: if the initial concentration of bacteria b 0 is not too large, coexistence depends weakly on b 0 , and very strongly on the initial algal concentration a 0 . For low enough bacterial concentrations, the smallest critical initial algal concentration for which coexistence will occur increases with λ. These features are reasonable considering that there is a diffusive delay in the metabolite exchange between reservoirs: if the delay is too long, auxotrophs will not be able to recover in the absence of a limiting nutrient.
The effect of the reservoir equilibration length η on the coexistence diagrams is more subtle. For small η, the crash-coexistence boundary sits above the membrane limit boundary, so toward  Figure 3: Transient dynamics of a bacterial population fed through a channel that allows metabolite diffusion from a remote carbon source. The diffusive exchange geometry controls the dynamics through the nondimensional channel length λ and reservoir equilibration length η. Model solutions predict that: (a) for fixed λ = 3, increasing η delays the time of peak bacterial growth and curtails growth due to a limited carbon-source flux; (b) for fixed η = 10, increasing λ significantly delays peak growth, with an impact on the maximum bacterial concentration attained. The delay as measured by τ max , the time of maximal growth rate, is proportional to λ 2 (inset). For all simulations, initial nondimensional bacterial and carbon concentration are b 0 = 5 × 10 −4 and c a (t = 0) = 10, other parameters are from table 3.
higher initial concentrations. This boundary is then pushed toward lower initial concentrations for intermediate values of η while still sitting above the membrane limit, before raising to higher initial values for high values of η. The general shape of the boundary is preserved for all η. To understand the nonmonotonic dependence of the boundary shift with η, we note η/λ is the reservoir/channel volume ratio. Thus, with λ fixed, changing η takes the populations through three regimes: i) the reservoir volume is small compared to that of the channel, η/λ 1; ii) the volumes are the same size, η/λ ∼ 1; iii) the channel volume is smaller than that of the reservoir, η/λ 1. In regime i), the equilibration time t eq = λη is small, but a large channel volume relative to the reservoirs dilutes any metabolite produced, making metabolites inaccessible to the microbial partner and preventing co-existence. In regime iii), the relative channel volume is small, but co-existence is impeded due to the long equilibration time t eq 1, which slows down significant metabolite exchanges between reservoirs. Finally, in regime ii), where reservoirs and channel have similar volume and t eq ∼ 1, mutualistic coexistence is favoured.
Aside from the co-existence or crash fixed points just discussed, we can use the model to analyse the transient dynamics leading to these equilibria. In particular, it is illuminating to evaluate the relaxation time taken for remote populations to reach the fixed points for a given initial microbial concentration in reservoirs assumed initially devoid of metabolites, as previously. Numerical solutions of the model equations show that this time varies as λ is increased across the co-existence/crash boundary for given η, as shown in Figure 5A. It is clear that the time to relax to the equilibrium rises sharply on either side of the critical λ at the boundary. This slow relaxation for λ values close to the bifurcation between extinction or co-existence is accompanied by oscillatory transients (see Figure  supplement 1). Similar considerations apply to the dependence of this time on the equilibration length η for a given λ, within that case there is the possibility of two boundaries between extinction and survival, see figure 5B.
Interestingly, the algal and bacterial concentration fixed points, a * , b * respectively, are independent of λ and η, see equations (2.22). Larger separation (increasing λ) or weaker diffusive coupling On the other hand, the response to an increase in η for fixed channel length λ = 3 is nonmonotonic. The coexistence initially contracts, then expands, and finally contracts again. The grey lines in both plots corresponds to the membrane limit for which λ → 0 and equilibration of metabolite concentrations between the two flasks is instantaneous. This provides the maximum possible concentration parameter space for mutualistic coexistence. The coexistence boundaries were determined by solving equations (2.14) and (2.15) numerically using the parameters in table 3. (c) Along the transect (dotted red line) in (a) corresponding to a conserved ratio of initial concentrations b 0 /a 0 = 20.0, the critical initial algal concentration a c 0 above which coexistence occurs is an increasing monotonic function of the length of tube λ. (d) Using the same transect in (b), the non-monotonic behavior of the critical algal concentration a c 0 with η is clearly revealed.
to the reservoirs (increasing η) delays the approach to equilibrium and reduces the extent of the mutualistic co-existence region. However, these geometric changes do not alter the microbial concentration fixed points, which have the same values as in the membrane limit. This equilibration is possible thanks to supply of metabolites (whose concentrations are also geometry-independent, see equations (2.22)) from the partner reservoir. A sufficiently large metabolite gradient across the channel is required to support the equilibrium metabolite and cell concentrations. Indeed, the model predicts an increase in the metabolite concentration at the production reservoir. For example, if the equilibrium concentration of carbon in the bacterial reservoir is c * b , then at the algal reservoir we predict c * where the function f can be obtained by comparison with equation (2.22). The same applies for the vitamins. This metabolite enrichment is an interesting prediction of the model. The concentration excess at the production reservoir is linear in both separation λ and equilibration length η: two parameters with which enrichment could be experimentally controlled. As an example, for the L. rostrata and M. loti mutualism using λ = 1.25 and η = 2 (all other parameters as before) our model predicts a sevenfold enrichment of vitamin B 12 in the bacterial reservoir compared to the algal side.

Discussion
The fundamental issue addressed here is how mutualistic dynamics depend on the geometry of diffusive exchanges. Phrased another way: under what conditions can one establish mutualistic interactions 'at a distance' ? To answer these questions, a mathematical model was solved to predict the behaviour of an algal-bacterial synthetic consortium, whose mutualistic dynamics in mixed culture have been experimentally characterised (Kazamia et al. 2012). Two key geometrical parameters control the diffusive exchange of metabolites: the separation λ (the nondimensional channel length) and the equilibration length η (the nondimensional ratio of growing volume to metabolite exchange area). Model solutions allow prediction of whether initial concentrations of algae and bacteria will result in mutualistic coexistence or population crashes for given values of λ and η. In particular, we obtain the boundary between regions exhibiting these two equilibria, whose shape is qualitatively unchanged when varying λ or η: it is approximately flat for a broad range of bacterial concentrations, falling very rapidly thereafter. The shape of the boundary can be intuitively explained by considering that an initially high concentration of one of the two species will produce a large initial amount of metabolite, which allows the partner species to grow and recover even from initially very low numbers. Thus coexistence will be achieved for abundant numbers of one or both species. When both species start at low concentrations, however, extinction occurs. The diffusive geometry parameters shift the boundary. We predict that increasing the channel length λ, at fixed equilibration length η and fixed bacterial concentration b 0 , increases the critical concentration of algae that will support co-existence with bacteria. On the other hand, for fixed λ and b 0 , the critical algal concentration varies nonmonotonically, falling and then rising again with increasing η. This behaviour can be explained in terms of the diffusive delays. The dependence on λ is intuitive: separating the partners further increases a diffusive delay, which we recall scales like λ 2 , so that more algae are required to support coexistence at a distance. The nonmonotonic behaviour with η is less obvious. It results from a dilution of metabolites in the volume of the channel for low values of η, requiring higher initial densities for successful coexistence, and from weak fluxes of metabolites into the homogenisation volume when η is large. With respect to these two extremes, coexistence is more easily achieved at intermediate values of η.
Contrary to the dynamics of a single microbial population sustained at a distance by a finite resource, for which the maximum achievable concentration depends strongly on the geometric coupling, our model predicts that mutualistic microbes at a distance can achieve as high a steady concentration as in a mixed environment. The diffusive geometry parameters only modify the transient dynamics and impose higher initial concentrations to avoid crashes. This result emphasises the need to consider that nutrient coupling between spatially segregated species over relatively large scales cannot be neglected when considering natural samples, for example sediment core communities.
There are many implications of this work for the microbial ecology of synthetic consortia. Indeed, a synthetic mutualism with one-dimensional diffusive exchanges could be engineered to test experimentally the predictions of the model. For example, the dynamics of two auxotrophic monocultures growing in separated flasks and linked by a channel allowing metabolite diffusion (filled with a hydrogel to prevent cross-contamination) could be quantified and compared with model predictions. In particular, it would be interesting to check if the metabolite accumulation at the production reservoir predicted by this work is experimentally verified. As mentioned, this would provide a method to enrich and chemically characterise trace metabolites, like vitamins. Preliminary experimental results on such connected flasks (see Materials and Methods) are consistent with this prediction and demonstrate the feasibility of diffusive coupling of mutualistic microbes in a macroscopic set-up. Alternatively, the model could be tested using cells in microfluidic chambers [e.g. as in the work of Kim et al. (2008)], although modifications would be necessary to account for stochastic effects associated with the small cell numbers in such systems (Khatri et al. 2012). Experiments with co-cultures exchanging metabolites across a membrane (Paul et al. 2013) could also explore the accuracy of the model in the membrane limit. As well as being tested, the model could be used to describe other synthetic consortia in which populations also interact diffusively across porous hydrogels (Kazamia et al. 2012, Harcombe et al. 2014 or microfluidic structures (Kim et al. 2008). It is straightforward to extend the model to account for two or three-dimensional diffusive exchanges appropriate to these systems.
The present model may also provide the foundation for a physical description of microbial networks, e.g. consortia for cooperative biosynthesis (Hom et al. 2015, Cavaliere et al. 2017 or microbial communities in soil, or spatially coupled biofilms (Liu et al. 2017). Indeed, as mentioned earlier, at the microbial scale, soil can be approximated as a physical network of growth chambers linked by channels (Pérez-Reche et al. 2012). In establishing the key geometric parameters that govern the most elementary unit in a network, namely two diffusively linked nodes (reservoirs), the present work provides a basis for describing population dynamics in a two-or three-dimensional network of coupled nodes (Figure 2). It is left to future work to take up the significant challenge of studying such networks, particularly when there is inhomogeneity in the diffusive couplings and stochasticity in the populations themselves. It seems likely that physical methods used to study lattice models in physics may prove useful in this problem. This view of microbial networks centring on the physics of diffusion could also help refine interaction matrix models of microbial communities and extend them beyond contact interactions (Mathiesen et al. 2011). Aside from the microbial networks mentioned above, the model may also be a relevant interpretative tool to understand the behaviour of structured environmental communities with diffusive exchanges, such as river biofilms (Battin et al. 2016) or sediment layers (Pagaling et al. 2014). Moreover, knowledge of the parameters for effective metabolite exchange between organisms that encounter one another is important to gain insight into how such communities initiate in the natural environment, and the drivers and constraints on the evolution of such mutualisms (Kazamia et al. 2016).

Diffusive reservoir equilibration (no microbes)
We consider here the purely physical equilibration between two diffusively connected reservoirs to reveal the interplay between the diffusive time and the equilibration time in such a system. This setup utilises the same geometry as in Fig. 2, with the reservoir atx = 0 having an initial concentration c 0 (t = 0) =c init of a chemical species, and the reservoir atx = L having an initial concentration c L (t = 0) = 0 of the same species. The chemical concentration along the tube is initially equal to zero, and has diffusivity D. Since our focus here is purely on the different physical timescales independent of biological processes, we choose a non-dimensionalisation scheme restricted to this section only that differs from the main body of the paper. Rescaling chemical concentrations by c init , lengths by L and time by L 2 /D, we obtain where we recognise the nondimensional parameter ζ = L/ , the ratio of tube length L to equilibration length = Γ/Σ. These equations are subject to initial conditions c 0 (0) = 1 , c L (0) = 0 , c(x, 0) = 0 and boundary conditions c 0 (t) = c(0, t) and c L (t) = c(1, t). Despite the fact that this is a linear PDE with apparently simple boundary conditions, the fact that it exists on a finite domain, and is coupled to the reservoir dynamics, makes it difficult to obtain an explicit analytical solution for general values of ζ.
Approximate solution for ζ 1 When ζ 1, the time evolution of the reservoir concentrations is much slower than the establishment of a concentration gradient in the tube. Thus, the diffusive dynamics within the tube reach a quasi-steady-state distribution between the two reservoir concentrations c 0 (t) and c L (t). In this approximation, the solution to the diffusion equation in the tube is the linear profile c(x, t) ≈ [c L (t) − c 0 (t)] x . Substituting this solution into the reservoir dynamics, and solving the resulting two ODEs yields (in dimensional units) (2.17) We thus deduce that in the limit ζ = L/ 1, the timescale of exchanges is purely dominated by the equilibration time τ eq = L /2D, as argued previously. The same time scale plays a role when the biological dynamics of growth and production are considered, as discussed in the main text. Finally, Laplace transforming the dynamical equations for the reservoir concentrations yields

General solution from Laplace transform
(2.20b) Combining the above we obtain explicit solutions for M (s) and N (s), thus entirely determining the solutions f 0 (s), f L (s) and f (s) to the problem in the Laplace space. In particular, for the concentration in the reservoir initially devoid of chemical, we obtain This solution in Laplace space is not easily inverted into an analytical expression for the evolution in time of c L (t) = L −1 (f L )|(t). In order to access its time evolution, we adapted a numerical inverse Laplace code in Python (Barbuto 2002) which implements the Zakian method (Halsted & Brown 1972, Abate & Whitt 2006. The numerical evaluation of c L (t), as a function of the characterisic nondimensional parameter ζ = L/ , is shown in figure 6. It reveals the typical nondimensional timescale of equilibration 1/2ζ, which in dimensional form becomes the previously discussed equilibration time τ eq = L /2D. At steady state, the concentration equilibrates between the two reservoirs and the tube at a final uniform value c f = 1/(2 + ζ). Finally, for ζ 1, the validity of the approximations of the concentration c L (t) as a saturating exponential in equation (2.17) is clearly demonstrated (Figure 6, right panel).

Mathematical model of remote mutualistic cross-feeding
Fixed points The dynamical system in equations (2.7) supports a trivial set of fixed points corresponding to a system with no cells (a = b = 0) and any residual concentrations of metabolites. The non-trivial fixed point is given by where we recall that σ j , with j = C or V, are the nondimensional metabolite production strengths; m i and κ i are the mortality/growth ratios and nondimensional uptake parameters for species i =A or B; λ = L/ b is the dimensionless channel length; η = / b is the dimensionless reservoir volume to area ratio or equilibration length; and = k a /k b and θ = D c /D v are the growth rate ratio and diffusivity ratio respectively. For the fixed point given by equations (2.22) to be physically relevant, the concentrations it describes must be positive. Therefore, the parameters must satisfy the following constraints: (2.23c) The first condition requires production strength to be strong enough to overcome cell mortality. This guarantees the existence of positive equilibrium algal and bacterial concentrations. The second and third conditions guarantee this positivity for carbon and vitamin concentrations, respectively. They require that microbial consumption be high enough to overcome production.

Membrane limit
The first natural limit of the model is that of zero channel length λ → 0, in which the reservoirs are in contact, but separated by a porous membrane. We call this the membrane limit because the membrane setup is as in membrane experiments (Paul et al. 2013), and we consider instantaneous equilibration of concentrations across the membrane as a good approximation. Fixed points for this limit are obtained trivially by letting λ → 0 in (2.22b)-(2.22c), which confirms that metabolite concentrations are equalised between reservoirs at steady state. We note that the membrane limit is identical to a mixed co-culture, where A and B grow mixed together in the same reservoir, except for the dilution effect associated with the segregation of the two species on either side of the membrane. The corresponding dynamical system for a mixed co-culture also admits a positive fixed point (a * , b * , c * , v * ) under the same conditions (2.23), with a * and b * given by (2.22a), c * = c * b from equation (2.22b) and v * = v * a from equation (2.22c). As mentioned earlier, such a co-culture model is fundamentally different from models considering mutualistic nutrient exchanges implicitly (Murray 1989, Yukalov et al. 2012, Grant et al. 2014, Holland & Deangelis 2010.

Remotely-fed monoculture
Another interesting limit is one in which a species in one of the reservoirs is replaced by a fixed concentration of metabolite. For example, we could have species B growing on C diffusing through the channel from the remote reservoir. In this limit, the model on the side of C reduces to passive diffusion from a source, which provides a useful control on the mutualistic dynamics, as mentioned in the results section. The mathematical model for such a remotely-fed monoculture is directly obtained from the remotely cross-feeding populations model (equations 2.7) by setting one microbial species and the metabolites it produces to zero.

Estimation of biological parameters
Monoculture experiments: Carrying capacities of M. loti and L. rostrata Liquid cultures of M. loti were grown for 3 days ( 33 • C, shaken at 240 rpm) in TY medium (tryptone 5 g L −1 , yeast extract 3 g L −1 , CaCl 2 · 2 H 2 O 0.875 g L −1 ) and washed in TP+ before serial dilution for counting of colony forming units. The post-wash concentration was estimated to be 5−10 × 10 8 cells mL −1 . Given the existing loss of cells during washing, we therefore allow the bacterial carrying capacity of our model K b to be in the range 5−50 × 10 8 cells mL −1 . Similarly, we estimated the carrying capacity of L. rostrata by growing these algae in TP+ with 100 ng L −1 of vitamin B 12 for 6 days to saturation (22 • C, shaken at 200 rpm, day/night cycle of 14h/10h), and plating them after washing in TP+ and serial dilution on TY agar plates for colony forming unit counting. We recorded saturation concentration ∼2 × 10 6 cells mL −1 , which, allowing for losses during cell washing, results in an accepted range of 1−10 × 10 6 cells mL −1 for the algal carrying capacity K a in our model.

Monoculture experiments: Death rate of M. loti
A pre-culture of M. loti in TY as above was washed in fresh TP+ and inoculated at a concentration b 0 = 3.2 × 10 8 cells mL −1 in 70 mL of TP+ without carbon source. Every two days, a 100 µL sample was taken to determine a live cell concentration through counting of colony forming units (CFUs) on TY agar. After a 2 days lag period, we measured an exponential decay of the bacterial population with death rate δ b ≈ 5 × 10 −2 h −1 over the next 6 days.

Co-culture experiments: Global fit of model parameters
The experiments whose outcomes were used to fit the model parameters utilised the following protocol. L. rostrata and M. loti were grown in TP+ medium at 25 • C on a 12h/12h day/night cycle, with 100 microeinsteins of light and shaking at 120 rpm. Bacterial concentrations were estimated with counts of CFUs on TY agar, and algal concentrations were obtained with a Coulter counter.
In some experiments, B 12 concentration was estimated with bioassays (Raux et al. 1996). Figure 7 shows the results for a set of six independent experiments (a-f) along with global fits to the model, corresponding to the values shown in table 3.

Mutualism at a distance: experimental proof of concept
To test experimentally the predictions of the mathematical model, we developed a system to culture mutualistic microbial species exchanging metabolites diffusively over a finite distance. Briefly, each of two 100 mL conical Erlenmeyer flasks was modified (Soham Scientific Ltd) to have a side arm (8 mm long, outside diameter 11 mm, inside diameter 9 mm) in which a small glass tube could be inserted (25 mm long, outside diameter 8.65 mm, inside diameter 7.45 mm). Sealing of the tube-flask junction was achieved by compression of O-rings on each side of a metal washer glued onto the glass tube (see figure 8a,b). The force of compression was established and maintained by mounting the flasks on custom sliding platforms (figure 8b,c). To prevent contamination, flasks were capped with silicon plugs (Hirschmann Silicosen type T-22) and aluminium foil, while the middle area of the flasks and tube assembly was also further covered with aluminium foil. The central glass tube connecting the inside of both flasks was filled with a polyacrylamide (PAM) gel (4% acrylamide w/v with a relative concentration of bis-acrylamide of 2.7%, filter-sterilised before pouring, BioRad). Once polymerised, the gels in their tubes were put in a bottle of sterile water and left to soak for 6 days to allow for any of the toxic non-polymerised monomer to diffuse out of the gel. We verified the very weakly hindered diffusion of B 12 through this gel by colorimetry, measuring a reduction of ∼ 10% of diffusivity with respect to B 12 diffusion in water, which validates the chosen gel pore size as allowing the diffusive transport of small metabolites. We also performed a test to check for cross-migration of the mutualistic species. Both flasks were filled with a rich bacterial medium for soil bacteria (TY), but only one side was inoculated with M. loti (see below for strain details). These bacteria reached a saturation density within a few days, but over a timescale of 2.5 months no bacteria were detected in the first flask, proving the PAM gel is not penetrable by bacteria (and by inference by the algae, which are larger). In such connected flasks, we inoculated one side with the B 12 -dependent green alga Lobomonas rostrata (SAG 45-1, wild type strain) and the other with the B 12 producing bacterium Mesorhizobium loti (MAFF 303099, wild type strain, original gift from Prof. Allan Downie, John Innes Centre, UK). Both inocula were diluted with TP+ medium (Kazamia et al. 2012) to the desired starting concentrations of microbes. The L. rostrata pre-culture was grown in TP+ with 100 ng L −1 of vitamin B 12 from colonies picked from a slant, while the M. loti pre-culture was grown in TY medium. Both pre-cultures were washed in fresh TP+ before inoculation in the assembly in order to remove any organic carbon and B 12 in the initial growth media. The initial concentrations of M. loti and L. rostrata were b 0 = 2.2 × 10 8 cells mL −1 and a 0 = 5.3 × 10 4 cells mL −1 , inferred from viable counts. To ensure culture sterility, flask assembly and inoculation were carried out in a laminar biosafety cabinet (PURAIR VLF 48). The connected flasks were mounted on a shaking platform (120rpm) within an incubator for 50 days, at 25 • C, with continuous illumination (80 µmol m −2 s −1 ). After this period, these assemblies were left in static incubation at 20 ± 2 • C and at ambient day/night light levels.

Viable counts and B 12 concentration measurements
Algal and bacterial populations were sampled 55 and 230 days after inoculation. No contamination (external or between species) was detected, and PCR screening was used to confirm species identity as Mesorhizobium loti bacteria and Lobomonas rostrata algae. This confirms the ability of the PAM gel to prevent cells from crossing, while allowing metabolites to be exchanged.
Viable counts revealed that the population of bacteria 55 days after inoculation was ∼ 10 3 smaller than the inoculum. At the same time point the algae had grown little: the cell concentration was only 1.3 times larger than the inoculum. After 230 days the bacteria had recovered, and the algae had grown significantly. At this time the algal concentration from two replicates was a = 7.8 ± 0.3 × 10 5 cells/cm 3 (where the uncertainty is the standard error in the mean), about 15 times the inoculation concentration and close to the carrying capacity they reach in well-mixed co-cultures (see table 1). While slight initial growth of the algae might be attributed to internal reserves of vitamin B 12 , it is difficult to account for growth 230 days after inoculation in the absence of the vitamin. Indeed, using bioassays (Raux et al. 1996) we measured a B 12 concentration of 24 ± 3 pg/ml in the medium on the side of the algae. On the side of the bacteria, we found 132 ± 7 pg/ml. This implies the existence of a concentration gradient across the tube between the two flasks. This is required for the supply of the B 12 to the algae, as predicted by the model (see equation 2.22c).  Figure 5. Example of oscillations of concentrations of cells and metabolites before convergence. Initial parameters are close to the boundary between survival and extinction (λ = 2, η = 3, a 0 = 2 × 10 −2 , b 0 = 3 × 10 −4 and no initial nutrients).