## Abstract

A mathematical model is presented for the growth of yeast that incorporates both dimorphic behaviour and nutrient diffusion. The budding patterns observed in the standard and pseudohyphal growth modes are represented by a bias in the direction of cell proliferation. A set of spatial indices is developed to quantify the morphology and compare the relative importance of the directional bias to nutrient concentration and diffusivity on colony shape. It is found that there are three different growth modes: uniform growth, diffusion-limited growth (DLG) and an intermediate region in which the bias determines the morphology. The dimorphic transition due to nutrient limitation is investigated by relating the directional bias to the nutrient concentration, and this is shown to replicate the behaviour observed *in vivo*. Comparisons are made with experimental data, from which it is found that the model captures many of the observed features. Both DLG and pseudohyphal growth are found to be capable of generating observed experimental morphologies.

## 1. Introduction

Yeasts are unicellular organisms belonging to the fungal kingdom that are typically between 4 and 40 micrometres in diameter. In contrast with moulds, which grow in elongated tubes of one or more cells called hyphae, yeasts typically grow as individual unconnected cells, sometimes referred to as the yeast form. Some species of yeast, designated dimorphic, can grow either in a colony of individual cells or in a linked chain of elongated cells, referred to as a pseudohypha, which resembles the tubular true hyphae of other fungi. The ale and baker's yeast *Saccharomyces cerevisiae* usually grows by budding in the yeast phase, but it can be induced to switch to the pseudohayphal phase under certain stress conditions [1]. Yeast species such as *Yarrowia lipolytica* and *Candida albicans* are naturally dimorphic and the switch from yeast-like to pseudohyphal growth is intrinsic to the latter species' activity as a human pathogen [2]. In addition, it has been argued that yeast cells are a useful model eukaryotic organism [3]. This diverse range of applications makes it of great interest to study the growth and, in particular, the morphogenesis of yeast cells.

Yeast colony morphology is controlled by the manner in which individual cells reproduce and the dimensions of the new cells produced. Typically, yeast cells reproduce asexually through mitosis, in which new cells are formed from protrusions, termed buds, that form on existing cells. In favourable growth conditions, such as a nutrient-rich environment, two patterns of bud-site selection are typically observed: the axial pattern and bipolar pattern [4], both of which are illustrated in figure 1. In the axial pattern, exhibited by haploid cells, the first bud forms near the bud scar, and subsequent buds form near previous budding sites, while diploid cells undergo bipolar budding, in which buds form at the opposite pole to that occupied by the previous bud.

Environmental conditions can trigger a switch in the growth patterns of both cell types, which alters the colony morphology. One such example is pseudohyphal growth, such as exhibited by the colony shown in figure 2. This growth pattern is characterized by three features: (I) a change to the distal-unipolar growth pattern; (II) the elongation of cells and (III) the adhesion of daughter cells to their mother. In the distal-unipolar pattern (figure 1) cells bud at the pole directly opposite the bud scar, while the cell elongation causes the axial aspect ratio to increase from approximately 1.5 up to 3.5 [5]. The chains of connected elongated cells that form during pseudohyphal growth are known as pseudohyphae as they resemble the single long hyphae cells of filamentous fungi. Pseudohyphal growth of *S. cerevisiae* is typically triggered by an external stress, such as the presence of fusel alcohols [6,7], which is thought to be a type of quorum sensing [1], heat [8] or low levels of mating pheromones [9]. Of particular interest here is pseudohyphal growth in response to nutrient limitation, which is a form of scavenging [1,5]. Through this response, non-motile yeast cells are able to search for nutrients away from the colony or invade a host organism [10].

Both experimental observations [11] and continuous coupled reaction–diffusion models [12,13] have shown that colonies of the motile bacterium *Bacillus subtilis* can produce growth patterns reminiscent of pseudohyphal yeast. Experimental studies [14–16] have shown that these patterns depend on both the nutrient concentration and the composition of the growth medium, and similar observations have been made for non-motile bacteria [15]. This suggests that bacterial colony morphology is influenced by diffusion-limited growth (DLG), which arises due to the interaction between cells and a second diffusing substance, such as nutrient or a byproduct of cell growth. It is known that yeast colony shape depends upon both the species of yeast and the limiting nutrient [17], while comparisons between *in vivo* data and a one-dimensional cellular automaton model coupled to a continuous nutrient model suggest that DLG alone does not control the colony morphology for all yeast species and nutrient types [18]. An off-lattice model has shown, however, that changes to the cell growth pattern are not able to alter the morphology without the inclusion of growth inhibition due to crowding [19].

Mathematical modelling clearly provides a promising avenue for identifying and quantifying the mechanisms that lead to filamentous growth. Discrete lattice-based models have been used extensively to understand pattern formation in cell colonies, such as the Eden model [20] and the diffusion-limited aggregation model [21], which represent growth in high and low nutrient concentrations, respectively. In the context of a growing colony, the latter is referred to as DLG. The interaction between growing cells and a diffusing species has been represented implicitly by probabilistic growth rules [22], through crowding effects [23], and using coupled cellular automata [24]. The effect of agar concentration and nutrient availability on yeast colony morphology has been investigated using coupled reaction–diffusion equations [25], while continuous models have also been used to analyse the stability of the colony front, albeit in the context of bacterial growth [13]. Using a hybrid model, comprising discrete non-motile cells on a lattice and a continuous nutrient field, it has been shown that DLG can be induced by increasing the nutrient consumption rate of the cells [26].

To understand the mechanisms that influence colony shape, we must be able to quantify the observed patterns. Bacterial colonies have been described by fractal dimension [11,14] and related scaling laws [12], while fingers growing from bacterial colonies have been quantified by the dimensionless width and amplitude [13]. Yeast invasiveness has been classified using the total colony volume and invasive colony volume [27,28], while a systematic analysis of 427 features identified six as important for classifying colony morphology, with the fractal dimension, average entropy texture within the colony and the area of the colony structure the most important features [29]. The size of the colony perimeter relative to that of its area and the coefficient of variation of the colony boundary as a function of angle have also been used to describe yeast colonies [25]. Much use has been made of normalized spatial counts [30,31] and pair correlation functions to describe spatial patterns of cell colonies [32–36]. In particular, normalized count functions and pair correlation functions have been used to define measures of cell distributions in the radial and angular directions for a growing colony of yeast cells [37].

While mathematical models have the potential to assist in the classification of yeast strains and mutants, this can only be realized if the simulated data can be compared with experimental data. Approximate Bayesian computation (ABC) allows the generation of probability distributions for parameters of interest [38,39], which provide more information than simple point estimates such as would be generated using a least-squares match. ABC algorithms require large numbers of simulations in order to generate the required distributions, which means that these techniques are only practical if the model solutions can be computed efficiently. This requires a careful balance between model accuracy and computational complexity.

Despite the substantial body of existing work, the models developed so far do not provide a sufficient understanding of the filamentous growth of yeast. In particular, no model incorporates all three key features (I, II and III from above) of pseudohyphal growth. Furthermore, there is little understanding of how pseudohyphal and DLG interact and, consequently, only a limited understanding of the conditions leading to filament-like patterns. We seek to address these questions using an efficient mathematical model for the dimorphic growth of yeast cells. Building upon the lattice-based discrete-time model developed by Matsuura [24], cells are represented by squares on a lattice, while a single species of nutrient is represented by particles that perform a random walk through the lattice. The changes in cell polarity and shape (properties I and II, respectively) that occur during pseudohyphal growth are modelled by modifying the original model of Matsuura [24] to include a bias in the direction of the cell growth, while cell–cell adhesion (property III) is enforced automatically by specifying that all cells are non-motile.

The model is used to examine the influence of nutrient diffusivity, nutrient concentration and the directional bias on the observed growth pattern. This was not possible using previous models, which did not incorporate both the diffusion of nutrient and the key features of pseudohyphal growth. The model thus offers a significant new insight into the relative impact on the morphology of each phenomenon, providing, for the first time, a theoretical understanding of the conditions under which each is important. The change in morphology that arises from pseudohyphal growth due to nutrient deprivation is investigated by increasing the bias as the surrounding nutrient concentration decreases. The morphology is measured using spatial statistics based upon those developed by Binder *et al.* [37], which have been shown to be suitable for use with experimental data. Using these statistics, we identify the dominant growth mechanisms in different regions of the parameter space. Importantly, the model developed here is efficient while still capturing important biological features, making it highly suitable for use with ABC algorithms and thus providing the first viable avenue for estimating physical parameters from experimental observations of pseudohyphal growth.

The remainder of this article is structured as follows. The mathematical model and statistics are described in §2 and the accuracy of the discrete nutrient model is discussed. In §3, the model is used to examine the influence of directional bias, nutrient concentration and nutrient diffusivity on the morphology, and the dominant growth mechanisms are identified. Some experimental comparisons are discussed in §4, while two extensions to the model are explored in §5. The main conclusions are summarized in §6.

## 2. Mathematical model

### 2.1. Model description

We adapt the single-layer diffusion model introduced by Matsuura [24], which is discrete in both time and space. Examining the experimental image in figure 2*c*, it is evident that a small section of the colony boundary, which grows in the radial direction, closely resembles rectilinear growth normal to the boundary. As inspection of this region is sufficient for determining the growth mode, we are thus able to determine the morphology by simulating only a small section of the boundary, which greatly reduces the number of computations required in the simulations. This approach avoids the difficulties of simulating a colony having a nearly circular boundary on a rectangular grid. The yeast cells are thus confined to a grid with *L*_{x} and *L*_{y} points in the *x*- and *y*-directions, respectively. The total area of the lattice is *A* = *L*_{x}*L*_{y}. Each point within the grid may be occupied by at most one yeast cell, with the number of cells at any time denoted by *ν* and the corresponding average cell density defined by *ρ* = *ν*/*A*. It is convenient to denote the initial number of cells *ν*_{0}, so that the number of new cells is *ν*_{g} = *ν* − *ν*_{0}. A single nutrient species is modelled by discrete packets containing mass *m* of nutrient that move on the same grid as the yeast cells. Each site of the lattice can hold any non-negative integral number of nutrient packets, regardless of whether that point is also occupied by a yeast cell. As an initial condition, we place a line of *L*_{x} cells across the row *y* = 1 of the grid, which matches the condition used by Matsuura [24] and represents a small portion of a colony boundary. We also specify the initial number of nutrient packets *n*_{0}, which are placed uniformly at random across the domain to give a total initial nutrient mass of *m*_{0} = *mn*_{0}. The initial average nutrient concentration is given by *c*_{0} = *m*_{0}/*A*, which is not necessarily uniform across the domain. The domain is taken to be periodic across the boundaries *x* = 1 and *x* = *L*_{x} but not in the *y*-direction. For a grid point (*x*, *y*) not on a boundary, we take the neighbouring cells to be the von Neumann neighbourhood consisting of (*x* + 1, *y*), (*x* − 1, *y*), (*x*, *y* + 1) and (*x*, *y* − 1). On a periodic boundary, the *y* co-ordinates are taken modulo *L*_{y}. On a non-periodic boundary, any point that lies outside the domain is ignored. At this juncture, we do not set values for the length, time and mass scales.

Given some state of the cells and nutrient, the state at the next time step is found by performing three operations: yeast cells may absorb nutrient; yeast cells may reproduce; and nutrient packets may move throughout the lattice. These mechanisms, illustrated in figure 3, operate as follows. Every yeast cell may absorb and store a nutrient packet with probability *p*_{a} provided that there is a nutrient packet on the same lattice point and the cell contains less than the maximum nutrient capacity *m*_{s}. Every yeast cell with at least mass *m*_{r} of nutrient stored attempts to reproduce with probability *p*_{r} by creating a new yeast cell in a neighbouring unoccupied grid location. As one neighbouring lattice site will always contain the mother cell, each cell can reproduce in three directions and hence can reproduce no more than three times. This limit reflects the crowding that occurs within a yeast colony. If the site selected for reproduction is occupied, then the reproduction event is aborted. When reproducing, the cell consumes *m*_{r} nutrient. Following Matsuura [24], we set *m*_{s} = 4; however, any choice for *m*_{s} is possible provided *m*_{s} ≥ *m*_{r}. The original model is modified so as to introduce a bias in the direction of cell reproduction that depends upon cell polarity. The proximal pole of the daughter cell is set to be the side adjacent to the mother cell, with the distal pole on the opposite side. The cell axis is defined to be the line joining these two poles. A reproducing cell creates a daughter in the site adjacent to its distal pole with probability *p*_{d}, and otherwise creates the daughter cell in either one of the two perpendicular sites with probability (1 − *p*_{d})/2, as illustrated in figure 3. Because no more than one cell can occupy a given site, a cell cannot create a daughter at the site adjacent to its proximal pole. The case corresponds to unbiased growth, while results in a bias towards growth in the direction of the cell axis, and a bias towards perpendicular growth. The exclusion rule means that each cell may reproduce no more than three times. At each time step, every nutrient packet walks randomly through the domain with probability *p*_{m}, and each selected packet moves *s* times at each iteration to one of the four adjacent sites with equal probability. We choose to vary *s* rather than *p*_{m} so as to allow a comparison with the original study by Matsuura [24]. Increasing *s* corresponds to increasing the diffusivity of the nutrient. The model parameters are summarized in table 1. Other features sometimes included in models of cell growth, such as metabolic states, cell death and subsequent lysis, and the production of cell waste products, are not included for reasons of simplicity.

### 2.2. Quantifying the spatial pattern

To describe the influence of the model parameters on the colony morphology, we require metrics that summarize the spatial pattern of the cells. Binder *et al.* [37] established that the morphology of disc-like yeast colonies could be quantified using scaled counts and a pair-correlation function computed from experimental data, which were summarized by three indices that measured: (i) the variation in the azimuthal direction; (ii) the variation in the radial direction; and (iii) the local aggregation of cells. Here, we develop analogues of the first two of these indices that are applicable to a rectangular geometry and provide a sufficient description of the pattern.

A colony of cells may be represented by a matrix *M* such that
so that the number of occupied sites is given by
with corresponding mean density *ρ* = *ν*/(*L*_{x}*L*_{y}). To describe the pattern, the cells are placed into the sets
with associated bin counts and .

To quantify the horizontal pattern, we first compute the standard deviation of *c*_{x}(*j*), given by
where *μ*_{x} is the mean of *c*_{x}(*j*). For non-exclusion processes, the standard deviation may be normalized by the value attained when all cells reside in a single bin [30]. While this has also been used as an upper bound for exclusion processes [31], if the number of cells is larger than the bin size *L*_{y}, then this upper bound can never be realized. If exclusion is considered, then a smaller upper bound on the standard deviation is
where *ν*_{y} = *ν* mod *L*_{y}. This corresponds to arranging the cells so as to create as many full bins as possible and placing any leftover cells in another bin. If the number of bins divides the number of cells exactly, then it is possible to have an equal number of cells in each bin, which corresponds to *σ*_{x} = 0; however, such an arrangement is not always possible in an exclusion process. In general, the minimum standard deviation is given by
where *ν*_{x} = *ν* mod *L*_{x}, which corresponds to arranging the cells as evenly as possible across the bins. Examples of these two distributions are shown in figure 4. When *L*_{x} divides *ν* exactly, this expression yields *σ*_{min} = 0. We then define the horizontal index
2.1This index provides a measure of the randomness in a given pattern and is an analogue of the azimuthal index of Binder *et al.* [37]. The index takes the value 0 for a pattern in which the cells are spread as evenly as possible across each value of *x*, and is unity when the cells are clustered in as few bins as possible.

As an analogue for the second index, we first define the total colony size *R* to be the maximum *y* index of the occupied cells, while the size *R*_{p} of the non-filamentous portion of the colony is defined to be the largest value of *y* at which *c*_{y} attains its maximum. The vertical index is then given by
2.2Larger values of *I*_{y} correspond to greater variation in the *y*-direction.

### 2.3. Accuracy of the discrete nutrient model

Compared to a continuum model, the discrete nutrient model employed here is advantageous as it facilitates the representation of the interaction between cells and nutrient packets, avoiding the need to couple the discrete cell model to a numerical solution method for a continuous diffusion model. In using this approximation, care must be taken to ensure the diffusion solution is resolved sufficiently so that the final cell distribution is independent of the number of nutrient packets. This is assessed using the spatial indices defined in §2.2.

In the original formulation of Matsuura [24], each packet in the model represents one unit of nutrient, which provides a relatively coarse approximation to the true nutrient diffusion. There is, however, no restriction on the nutrient content of each packet. If the total starting nutrient mass *m*_{0} is spread evenly across *n*_{0} nutrient packets, the mass of each packet is
which may take on rational values. The original model of Matsuura [24] corresponds to the case *n*_{0} = *m*_{0} and, for any fixed value of *m*_{0}, increasing *n*_{0} provides a better approximation to the continuum limit.

To investigate the accuracy of the discrete model, we first consider a representative problem with nutrient diffusion only and compare the solutions from the discrete model to the continuous analogue. For the discrete problem, we consider a domain with dimensions *L*_{x} = *L*_{y} = 100 initially populated randomly with nutrient at lattice sites with coordinates satisfying 46 ≤ *x* ≤ 55, with total mass *m*_{0} = 10^{5}. The number of steps is nominally set to *s* = 1. The solution at each iteration is averaged over all values of *y* to provide a one-dimensional representation of the diffusion.

The corresponding continuum problem has domain half-width *L* = 49.5, while the populated region has half-width *a* = 5. When *s* = 1 it is known that the nutrient concentration at position *x* and time *t*, denoted by *C*(*x*, *t*), is governed by the linear diffusion equation [40]. When *s* > 1, each nutrient packet can undergo multiple motility events at a given time step, which simply corresponds to a change in the diffusion timescale. Taking this into account by rescaling time, the continuum model for general *s* is governed by
where *D*_{n} = *sp*_{m}/4 is the coefficient of diffusion.

The packet concentrations from the discrete model calculated using different initial packet numbers *n*_{0} are plotted in figure 5*a*, along with the initial condition and series solution to the continuum model. As *n*_{0} increases, the discrete solutions approach the continuum solution and, by *n*_{0} = 10^{6}, the discrete and continuum solutions are indistinguishable at the plotted resolution.

While this analysis suggests that a large number of nutrient packets are required in order to find an accurate solution to the diffusion model, we are, ultimately, only interested in determining the values of the indices *I*_{x} and *I*_{y}. To investigate the dependence of these on *n*_{0}, we consider the full model with both cell growth and nutrient diffusion for unbiased growth with *s* = 2 and *m*_{0} = 4 × 10^{4} on a lattice with *L*_{x} = *L*_{y} = 100, corresponding to the initial concentration *c*_{0} = 4. The model is run up to *ρ* = 0.5 using values of *n*_{0} between 40 and 4 × 10^{6}. For each value of *n*_{0} we compute 50 realisations and calculate the averages of the associated indices. To facilitate a comparison between different values of *n*_{0}, the cells are allowed to absorb multiple nutrient packets at each iteration up to a maximum total mass of 1. The value of each index is plotted against *n*_{0} in figure 5*b*, from which we observe that both indices are approximately constant for *n*_{0} ≥ *m*_{0}. The mean values of both indices at *n*_{0} = *m*_{0} have relative errors of no more than 0.06 compared to the values found with the largest number of packets *n*_{0} = 4 × 10^{6}. This indicates that increasing the accuracy of the discrete nutrient model beyond *n*_{0} = *m*_{0} does not result in a significant change in the spatial indices.

While results for only one set of parameters has been shown here, we note that the same outcome is found when using other values for *s* and *c*_{0}. We thus conclude that accurate values of *I*_{x} and *I*_{y} may be obtained for relatively small values of *n*_{0}. In particular, setting *n*_{0} = *m*_{0} provides a good estimate of the indices and this choice will be used for the remainder of this study.

## 3. Directional bias and morphology

### 3.1. Unbiased growth

As a reference point, we first consider colony growth in the absence of directional bias, which corresponds to setting and represents bipolar growth. Following previous experimental [16] and mathematical studies [24,25], we investigate how both the initial nutrient level and the nutrient diffusivity affect the morphology. We take *L*_{x} = *L*_{y} = 100 and consider initial nutrient concentrations *c*_{0} between 1 and 7, and nutrient step numbers *s* varying from 1 to 37, which are chosen to best illustrate the observed behaviour. For convenience, we refer to each solution by the corresponding ordered pair (*s*, *c*_{0}). In each case, the solution is computed until *ρ* = 0.5, or until all the nutrient is consumed, so as to provide a comparison of the morphology between each set of parameters. While there is enough nutrient for the colony to reach *ρ* = 0.5 for simulations with *c*_{0} > 1.47, in practice some nutrient may be stored within cells that are unable to reproduce, which means that the simulation cannot reach the threshold density. To illustrate the general behaviour for each set of parameters (*s*, *c*_{0}), we compute 50 realizations.

Representative morphologies for each pair (*s*, *c*_{0}) are shown in the phase diagram in figure 6. For *s* = 1 and *c*_{0} ≤ 2, the colony grows into thin branches reminiscent of DLG growth. Significantly, this suggests that the colony can undergo filamentous growth without any bias in the direction of cell growth. This behaviour was not shown in the corresponding diagram given by Matsuura [24, fig. 6, p. 318] for the multilayer version of the model. For 5 ≤ *s* ≤ 17 and *c*_{0} ≤ 5, the colony grows in thicker branches. If either *s* or *c*_{0} becomes sufficiently large, then the colony grows uniformly, irrespective of the value of the other parameter.

The changes in morphology are quantified by the spatial indices *I*_{x} and *I*_{y}. To illustrate the general behaviour of the colony growth we compute the mean value of each index over 50 realizations. The mean values of the indices over these realizations are denoted by and , respectively, and are plotted in figure 7. The indices show that the greatest filamentous growth, occurs at low values of *s* and *c*_{0}, and that this decreases as either *s* or *c*_{0} increases. While the indices all share this general trend, the rate of change of each with respect to *s* and *c*_{0} is not the same, which indicates that each index provides unique information about the spatial pattern. By providing a quantitative measure of the morphology, the plots in figure 7 represent an improvement on previous phase diagrams, which classify morphology qualitatively, based on observations, rather than measurements [15,16].

### 3.2. Biased growth

To illustrate the effect of directional bias on the morphology, we consider two extreme cases: *p*_{d} = 0.2 and *p*_{d} = 0.8. In the former case, daughter cells grow perpendicular to the cell axis, which runs between the proximal and distal poles, 80% of the time, similar to the axial growth pattern; whereas, in the latter, daughter cells grow from the distal pole 80% of the time, which resembles the distal-unipolar growth pattern. Solutions are computed for the same values of *s* and *c*_{0} as in the unbiased case.

Typical colony shapes for *p*_{d} = 0.2 are plotted in figure 8*a*. Notably, simulations with *s* ≤ 13 and *c*_{0} ≤ 4 produce filamentous patterns despite a preference for perpendicular growth. This behaviour is a form of DLG, similar to the behaviour observed in non-motile bacteria at low-nutrient concentrations [16]. Representative colonies for *p*_{d} = 0.8 are shown in figure 8*b*. The large value of *p*_{d} results in straight filaments that are typically longer than those for unbiased growth.

The changes in behaviour due to varying *p*_{d} are quantified by the mean values of the indices *I*_{x} and *I*_{y}, plotted in figure 9, computed using 50 realizations for each pair (*s*, *c*_{0}). By inspecting the indices, the changes in morphology observed between the two values of *p*_{d} can be placed broadly into three categories. When both *s* and *c*_{0} are small, the colonies display large index values regardless of *p*_{d}, such as at (*s*, *c*_{0}) = (1, 1). Similarly, large values of both *s* and *c*_{0} generate colonies with low indices, as seen for (*s*, *c*_{0}) = (37, 7). Within the intermediate region the growth pattern depends strongly on the bias *p*_{d}. For example, colonies with (*s*, *c*_{0}) = (21, 3), show a significant increase in both indices between *p*_{d} = 0.2 and 0.8. Thus, morphologies in this region may be altered by adjusting the directional bias. This matches the dimorphic behaviour of real-world yeast colonies, which may exist within this parameter regime.

These three growth regimes may be delineated using one of the spatial indices. We first define to be the value of *I*_{x} for distal-unipolar (biased) growth divided by the average value computed for bipolar (unbiased) growth at the same parameter values (*s*, *c*_{0}). The three growth regimes are then defined by the following criteria:
3.1a
3.1b
3.1cThese three regions, plotted in figure 10, partition the parameter space with respect to *s* and *c*_{0} based upon how each parameter pair is affected by the change to distal-unipolar growth. Within the uniform and DLG regions the morphology is relatively insensitive to the cell growth pattern, remaining uniform or filamentous, respectively, regardless of the growth pattern. The intermediate zone consists of a region between the DLG and uniform zones in which the morphology varies significantly with the growth pattern.

As discussed above, all simulations are computed until the cell density *ρ* reaches 0.5 or all nutrient is consumed. All of the colonies that do not achieve the target density lie in the region *c*_{0} ≤ 3. In each case, the mean density increases with the bias *p*_{d}, meaning that colonies with a greater bias produce a larger number of cells. This illustrates that the growth mode of the cells has an impact on the survival of the colony, and thus the ability to alter this growth mode provides an evolutionary advantage.

## 4. Experimental case studies

### 4.1. Methodology

We consider two experimental case studies: the onset of filamentous growth and very mature colony growth. For each study, we seek to infer the dominant growth mode in the experimental data by comparing the morphology with the behaviour predicted by the model. We emphasize that these comparisons are to illustrate general behaviour only and rigorous comparisons of model results with experimental data are left for future work.

A yeast colony was seeded from a single cell of *Σ*1278b (diploid, prototrophic) on 55 mm plates with yeast nitrogen base (YNB) agar, at a depth of 2 mm (Difco^{TM} YNB w/o amino acids and ammonium sulfate, cat. no. 233520, with 2% glucose, 50 μM ammonium sulfate and 2% twice washed bacto agar). Plates were incubated at 30°C. Each image is a typical representation of three biological replicates.

### 4.2. Onset of filamentous growth

As an example of the onset of filamentous growth, we consider the experimental image from figure 2*c*, which was taken 13 days after inoculation with a single cell, at which time the colony had a diameter of approximately 3.2 mm. To identify the growth mechanism that generated this pattern, we compare this image to simulations of both unbiased growth () and biased growth (*p*_{d} = 0.8). To achieve sufficient resolution in the simulations, we set *L*_{x} = 300 and *L*_{y} = 100, and 50 realizations were computed over the same range of parameters as in §3. As the simulation represents only a small section of the boundary, it is not appropriate to make a comparison using the value of *I*_{y}, which must be computed from an entire colony. The data and simulations are hence compared using *I*_{x} alone, with the experimental image yielding *I*_{x} = 0.266. For each value of the bias *p*_{d}, we rank the simulations by how closely *I*_{x} matches the experimental value, select the top 100 cases and identify to which of the three regions shown in figure 10 they correspond.

With unbiased growth, the best match is provided by (*s*, *c*_{0}) = (1, 4) in the uniform region, for which *I*_{x} differs by less than 1.2 × 10^{−5}, while the top 100 matches all lie within 8.5 × 10^{−3} of the experimental value. We find that 81% of these matches come from the DLG region, 7% are from the intermediate region and 12% are from the uniform region. For the biased example, the best match is given by (*s*, *c*_{0}) = (29, 3), which lies in the intermediate region, and differs from the experimental value by only 8.6 × 10^{−5}. In contrast with the unbiased example, the top 100 matches are more diverse, comprising 18% DLG growth, 46% intermediate growth and 36% uniform growth.

While the composition of the top 100 results differs between the unbiased and biased examples, each of these is able to produce similar morphologies. To illustrate this, we compare the best match from the DLG region computed using unbiased growth, which was the fifth-best match overall for unbiased growth, and the best match from the intermediate region computed using biased growth, which was also the best match overall for unbiased growth. These simulations had parameters (*s*, *c*_{0}) = (13, 3) and (*s*, *c*_{0}) = (29, 3), respectively, and are plotted in figure 11. Both simulated images show similar behaviour and have values of *I*_{x} within 1.5 × 10^{−4} of the experimental image. Consequently, this model suggests that it is impossible to tell DLG from pseudohyphal growth by examining images using *I*_{x} alone.

### 4.3. Very mature colony growth

As a second case study, we consider the morphology of a colony after 58 days of growth under the same conditions as previously, shown in figure 12*a*, which represents a very mature yeast colony. As a model comparison we use a domain with *L*_{x} = *L*_{y} = 200 with seed cells placed at every grid point within 60 units of the domain centre, which represents a region of uniform growth. This initial condition reduces the number of computations required and avoids the issues inherent in simulating colonies with circular boundaries on a rectangular grid. We consider unbiased growth () and consider a range of parameters (*s*, *c*_{0}), in each case computing the simulation until *ρ* = 0.5 or all nutrient is consumed. This value is chosen so as to provide a similar level of growth to the experimental image. We aim only to compare the growth patterns produced by the experimental colony and the simulations, rather than make a precise match.

To compare the model simulations with the experimental image, we use two of the spatial indices developed by Binder *et al.* [37] for disc-like colonies that motivated the indices *I*_{x} and *I*_{y}; respectively, these are *I*_{θ}, which measures angular variation, and *I*_{r}, which measures radial variation. To match the definition of *I*_{x} we here define *I*_{θ} to be the square root of the version used by Binder *et al.* [37]. For both indices, the data are grouped into 50 bins, which has been found to provide sufficient resolution, with radial bins used to compute *I*_{r} and angular bins used for *I*_{θ}.

The experimental image in figure 12*a* has indices *I*_{r} = 0.29 and *I*_{θ} = 0.113. Computing a single realization of the model using (*s*, *c*_{0}) = (1, 2) produces a colony with *I*_{r} = 0.27 and *I*_{θ} = 0.119, while a single realization with (*s*, *c*_{0}) = (1, 3) produces a colony with *I*_{r} = 0.23 and *I*_{θ} = 0.0936, both of which are plotted in figure 12. Each simulation closely resembles the experimental image qualitatively, while the index values show that these also provide a good quantitative match. Both sets of parameters lie within the DLG region, while simulations outside of this region produce smaller indices that do not provide a good match to the experimental image. Further, simulations with biased (pseudohyphal) growth are able to produce similar patterns. This suggests that, when measured by *I*_{r} and *I*_{θ}, it is not possible to distinguish the filamentous growth of the colony in figure 12*a* from a morphology produced by DLG alone.

## 5. Extensions

### 5.1. Ongoing nutrient consumption

As in the original formulation of Matsuura [24], it has been assumed that the cells require a specified amount of stored nutrient so as to reproduce but do not need to consume nutrient to survive; however, in practice, cells require ongoing nourishment for maintenance [41]. The consumption of nutrient may be modelled by removing the limit *m*_{s} on the total amount of nutrient each cell can store. Typically colony shapes with this restriction removed are shown in figure 13, while the mean indices from 50 simulations for a range of values for *s* and *c*_{0} are plotted in figure 14. This change in *m*_{s} causes increases in and , plotted in figure 13, compared to the values obtained for limited nutrient consumption from figure 6. This occurs because removing the storage limit reduces the amount of nutrient available over the course of the simulation, which leads to DLG and hence greater non-uniform behaviour across almost all parameters.

By computing simulations with ongoing nutrient consumption and *p*_{d} = 0.8, we are able to identify the three growth regions (3.1), which are shown in figure 15. Including ongoing nutrient consumption causes a shift in the location of the three growth regions identified in §3, and thus must be considered when predicting experimental results based on known physical parameters. We emphasize, however, that ongoing nutrient consumption does not affect the existence of the three growth regions from §3 but only alters the values of *s* and *c*_{0} at which they occur. While the model could be further extended to include the degradation of stored nutrient by the cells, this would again only alter the location of the growth regions and hence is not considered here. Future models that include ongoing nutrient consumption should incorporate nutrient degradation in order to provide a more realistic representation of this process.

### 5.2. Bias due to nutrient concentration

The examples considered thus far have demonstrated that increasing the directional bias causes some colonies to change from uniform to filamentous growth. In all of these examples, the bias has been held constant throughout the simulation; however, in practice the transition to pseudohyphal growth of a yeast colony is typically attributed to nutrient deprivation, which can occur either when the medium contains little nutrient or when nutrient cannot diffuse with sufficient speed to feed the colony. This transition may be modelled by taking the initial bias to be , corresponding to uniform growth, and increasing this value as the local nutrient concentration falls. For the purposes of this study, we define the local concentration *c* at a given cell to be the average concentration within the lattice points that can be reached in at most two steps, as illustrated in figure 16.

We here consider two forms for this change. The first is a linear transition from unbiased growth () at nutrient concentrations above some critical concentration *c*_{t} to biased growth (*p*_{d} = 1) at local concentration *c* = 0. This corresponds to setting
In the second form, the growth is uniform () for local concentrations *c* > *c*_{t} and jumps to a bias of 0.8 for 0 < *c* < *c*_{t}, which is described by
To facilitate a general investigation of this transition, the value of *c*_{t} is chosen to be the minimum amount of nutrient required to reach *ρ* = 0.5 with *L*_{x} = *L*_{y} = 100 and *ν*_{0} = 100, so that
In practice, the value of *c*_{t} depends on the strain of yeast used and must be determined experimentally.

Typical colony shapes for the two bias rules are shown in figure 17 for a range of initial nutrient concentrations *c*_{0} and nutrient steps *s*, with cells coloured by the value of *p*_{d} of their mother cell at the time of birth. Comparing the two cases, it is clear that the linear and step cases result in noticeably different biases at each parameter pair (*s*, *c*_{0}), as indicated by the differences in colour between the two figures. The mean indices and are computed from 50 realizations of each pair of parameters and are plotted in figure 18. Both indices indicate that similar patterns are produced by both bias rules, with the most significant difference being a larger value of at (5, 2) for the step rule compared to the linear rule. This suggests that both the linear and step bias rules produce similar colony morphologies, despite the differences in the bias *p*_{d} throughout the colonies.

Both of the bias rules cause increases in the indices at intermediate values of *s* and *c*_{0} compared to the values for uniform unbiased growth plotted in figure 7. This is in agreement with the results from §3, where it was found that the growth pattern is strongly dependent upon the bias at intermediate values of *s* and *c*_{0}. The increase in the indices also agrees with the experimental observation that nutrient deprivation leads to pseudohyphal growth.

## 6. Conclusion

We have developed a two-dimensional agent-based model for the growth of yeast cells that incorporates both dimorphic growth behaviour and nutrient diffusion. In this model, non-motile yeast cells occupy sites on a lattice and absorb discrete nutrient packets that walk randomly throughout the same lattice. The growth pattern is specified by adjusting the bias in the direction of cell proliferation. Improving on previous work, this model accounts for all three key features of pseudohyphal growth: (I) a change to the distal-unipolar growth pattern; (II) the elongation of cells and (III) the adhesion of daughter cells to their mother. The model developed captures the important biological features of pseudohyphal growth while remaining computationally efficient, and hence is highly suitable for use with statistical inference methods, such as ABC algorithms, to compute parameter estimates from experimental data. The spatial patterns that develop were quantified using two indices *I*_{x} and *I*_{y}, which are related to similar indices that have been used by Binder *et al.* [37] to quantify experimental images of yeast colonies. Larger values of these indices indicate stronger filamentous patterns. Simple rules have been used to link the bias to the nutrient level, which replicates the real-world behaviour of yeast colonies.

In general, increasing the bias causes an increase in both indices, indicating that the bias is correctly capturing the transition to pseudohyphal growth. By varying the initial nutrient concentration and diffusivity, and comparing the average value of *I*_{x} computed for a number of realizations, the parameter space may be divided into three regions: the uniform region, in which the colony front is smooth; the DLG region, in which the front is always rough; and the intermediate region, in which the morphology depends strongly upon the bias level. If the transition to pseudohyphal growth is solely due to budding pattern and cell length, then real-world colonies probably exist in this region. Outside of the intermediate region the morphology is largely determined by the nutrient concentration and diffusivity.

Comparing simulations of unbiased and biased growth showed that both DLG and pseudohyphal growth can produce similar patterns, and both mechanisms can provide a good match to experimental images when measured using the spatial indices employed here. This suggests that both mechanisms may contribute to the morphology and that it is not possible to distinguish DLG from pseudohyphal growth using the indices alone without knowing the precise growth conditions. The experimental results of Chen *et al.* [25] show that yeast colony morphology changes as the agar concentration varies between 1.5% and 6% with glucose (nutrient) levels between 0.5% and 2%. The diffusion of sodium and caesium has been measured in gels of up to 4% agar at 10°C and 25°C, from which it was found that the diffusion coefficient could be approximated by *D*_{n} = *D*_{0}(1 − 0.023*ω*), where *D*_{0} is the diffusion coefficient in water and *ω* is the weight percentage of agar [42]. This suggests that the diffusivity varies little with agar concentration. There is thus a need for a greater understanding of the role that nutrient diffusivity plays in yeast colony morphology; in particular, how this varies with agar concentration, and how this interacts with pseudohyphal growth patterns. The model developed here provides an avenue for exploring these questions, which was not possible using previous models. Furthermore, there is scope to use this model to understand the growth of other microbial colonies and to compare the dominant mechanisms leading to non-uniform morphologies, providing a complete theory of these processes.

The degree of pseudohyphal growth in the model may be adjusted by varying the directional bias parameter. By tuning this parameter using experimental data and, for example, ABC algorithms, the model has the potential to facilitate the identification of strain-specific growth characteristics, paving the way for the development of strain-classification schemes. The identification of strain-specific features is also vital to the analysis of gene-deletion assays. Furthermore, the greater understanding of yeast cell growth provided by this model allows the cell behaviour to be explored and optimized, either to encourage growth, such as is desired for food production, or to prevent it, such as in the case of yeast colony growth on medical equipment.

While the model considered here captured the general behaviour of yeast colony growth, all of the patterns were constrained to a lattice. This limited the number of daughter cells to three and constrained the patterns that could be realized. In practice, cells can produce approximately 24 daughters, and approximately seven when confined within a colony, and bud in a variety of directions [43]. In addition, the budding pattern of pseudohyphal growth and the corresponding elongation of cells have both been represented by an increase in the directional bias. These factors are likely to have an influence on the colony morphology. Future work will consider off-lattice models with fewer restrictions on cell growth and that incorporate changes in cell shape explicitly. Consideration must also be given to the effect on the nutrient level of lysis products from dead cells.

Both DLG and pseudohyphal growth appear similar when measured by *I*_{x}; however, it is possible that a more sophisticated measure, including a combination of different metrics, could distinguish these growth modes. Future work will thus continue to investigate and develop suitable spatial measures, along with techniques for categorizing colony morphology using these measures.

## Data accessibility

The processed experimental image analysed in figure 11 is included in the electronic supplementary material. The code used to generate the simulations is available in the online electronic supplementary material.

## Authors' contributions

B.J.B., J.M.G., J.F.S., S.G.O. and V.J. conceived and designed the experiments. J.M.G. performed the experiments and collected the data. H.T. and B.J.B. undertook the mathematical modelling and data analysis. H.T., B.J.B., J.M.G., J.F.S., S.G.O. and V.J. wrote and edited the paper. All authors gave final approval for publication.

## Competing interests

We declare that we have no competing interests.

## Funding

H.T. and B.J.B. were supported by Australian Research Council Discovery Project DP160102644, awarded to B.J.B. and S.G.O. V.J., J.M.G. and J.F.S. were supported by Australian Research Council Discovery Project DP130103547 awarded to V.J. and S.G.O.

## Acknowledgements

The *Σ*1278b yeast strain used in the experiments was supplied by the Fink Laboratory at the Whitehead Institute for Biomedical Research.

## Footnotes

Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3876076.

- Received April 30, 2017.
- Accepted August 31, 2017.

- © 2017 The Author(s)

Published by the Royal Society. All rights reserved.