Host control and nutrient trading in a photosynthetic symbiosis

Photosymbiosis is one of the most important evolutionary trajectories, resulting in the chloroplast and the subsequent development of all complex photosynthetic organisms. The ciliate Paramecium bursaria and the alga Chlorella have a well established and well studied light dependent endosymbiotic relationship. Despite its prominence there remain many unanswered questions regarding the exact mechanisms of the photosymbiosis. Of particular interest is how a host maintains and manages its symbiont load in response to the allocation of nutrients between itself and its symbionts. Here we construct a detailed mathematical model, parameterised from the literature, that explicitly incorporates nutrient trading within a deterministic model of both partners. The model demonstrates how the symbiotic relationship can manifest as parasitism of the host by the symbionts, mutualism, wherein both partners benefit, or exploitation of the symbionts by the hosts. We show that the precise nature of the photosymbiosis is determined by both environmental conditions (how much light is available for photosynthesis) and the level of control a host has over its symbiont load. Our model provides a framework within which it is possible to pose detailed questions regarding the evolutionary behaviour of this important example of an established light dependent endosymbiosis; we focus on one question in particular, namely the evolution of host control, and show using an adaptive dynamics approach that a moderate level of host control may evolve provided the associated costs are not prohibitive.


Abstract
Photosymbiosis is one of the most important evolutionary trajectories, resulting in the chloroplast and the subsequent development of all complex photosynthetic organisms. The ciliate Paramecium bursaria and the alga Chlorella have a well established and well studied light dependent endosymbiotic relationship. Despite its prominence there remain many unanswered questions regarding the exact mechanisms of the photosymbiosis. Of particular interest is how a host maintains and manages its symbiont load in response to the allocation of nutrients between itself and its symbionts. Here we construct a detailed mathematical model, parameterised from the literature, that explicitly incorporates nutrient trading within a deterministic model of both partners. The model demonstrates how the symbiotic relationship can manifest as parasitism of the host by the symbionts, mutualism, wherein both partners benefit, or exploitation of the symbionts by the hosts. We show that the precise nature of the photosymbiosis is determined by both environmental conditions (how much light is available for photosynthesis) and the level of control a host has over its symbiont load. Our model provides a framework within which it is possible to pose detailed questions regarding the evolutionary behaviour of this important example of an established light dependent endosymbiosis; we focus on one question in particular,

Introduction
Endosymbiotic relationships are widespread in nature and play key roles in the functioning of many ecosystems [1,2,3,4,5,6,7]. Different symbioses have evolved many times throughout history; of particular note is the evolution of cellular organelles such as chloroplasts from a cyanobacteria-eukaryote symbiosis [8]. One well-known example of endosymbiosis is the relationship between the ciliate Paramecium bursaria and the alga Chlorella [9]. These organisms [9,10] and their close relatives [1,11,12] have been the focus of much study in both symbiotic and free-living contexts, and therefore provide an excellent model system for the study of alga-protist endosymbioses.
The primary benefit of an endosymbiotic relationship between a heterotophic host and a photosynthetic symbiont (photosymbiosis) is thought to be nutrition [4,5,9]. The host obtains nutrients from its environment via phagotrophy-the engulfing of cells or particles and subsequent digestion within a vacuole. Free-living algae are also ingested in this manner, but not all digested; rather, some resist digestion long enough for a section of the digestive vacuole membrane to 'pinch off' and form a new, distinct vacuole. Known as the perialgal vacuole, this provides an alga with protection from digestion [13,14], allowing it to take up residence within the host ciliate and carry out the usual unicellular life cycle of growth and cytokinesis (division). Such symbiotic algae are now dependent on their host for nutrients which are unobtainable via photosynthesis (in particular, nitrogen). In return for these nutrients, the symbiont releases a portion of its photosynthate into the host cytoplasm, resulting in a net gain of organic carbon for the host [15,16,17].
The consequence of this nutrient exchange is that the photosymbiosis exists on a context-dependent continuum whereby the nature of the interaction depends on the light level [36]. In low light, the correspondingly low photosynthetic output of the symbionts results in a net loss of nutrition for the host-this is effectively parasitism. As light increases, the increase in pho-tosynthesis results in symbionts providing a net nutritional benefit to their host, yielding a mutualistic relationship.
A key step in the establishment of a permanent symbiotic relationship is the maintenance of a stable symbiont population [18]. Clearly, if the host population grows more rapidly than the symbiont population, successive generations of host cells will become increasingly diluted until a completely aposymbiotic state is reached. This can occur when P. bursaria are grown in the dark [9]. Conversely, if the symbiont population is the faster growing of the two, it will increase to the point of saturation, with severe consequences for the host-the symbionts have become parasites. Hence, if a host is to maintain a stable symbiont population it must carefully balance the gain and loss of symbionts so the two populations increase at approximately equal rates, either through controlling the rate of intake of new symbionts and the rate of removal (through expulsion or digestion), or by synchronising the cellular division cycles of the organisms. Both of these mechanisms could potentially lead to conflict between host and symbiont and the need for greater control by the dominant partner-presumed to be the host-to maintain the symbiosis.
Alga-protist endosymbioses have been addressed only briefly in the mathematical literature; see [19] for a review. There has been much emphasis on potential mechanisms for cell-cycle synchronisation in Chlorella-Hydra symbioses [20,21] and Chlorella-ciliate symbioses [22], while others focus on the role of nutrient trading [23,24] in more general photosymbioses. Conditions determining the evolution of an obligate endosymbiosis were investigated in [25,26,27]. A relevant recent paper modelled the P. bursaria-Chlorella symbiosis, showing how the combination of dynamic nutrient trading and differences in growth rates between partners yield a steady symbiotic population, but neglected to incorporate the potentially significant effects of symbiont intake and removal [28]. To the best of our knowledge, the present work represents the first attempt to provide a comprehensive ecological model encapsulating the myriad facets of symbiosis across a range of environmental conditions.
In this article we develop a model to illustrate the mechanistic basis for a photosymbiotic relationship and the configuration of the resultant mixotrophic holobiont. We describe host-symbiont interactions by a deterministic system of ordinary differential equations, incorporating the vertical transmission of symbionts via host cytokinesis and the horizontal transmission of symbionts via ingestion from, and egestion into, the environment. The above discus-sion on host-symbiont cell-cycle synchronisation forms the basis of a key assumption; namely that on the timescale of our model, host and symbiont cytokinesis is almost concurrent, in that daughter cells have a symbiont load equal to that of their mother cell. This has been directly observed in P. bursaria [29], and is in contrast to asynchronous cell cycles, for example, in which daughter cells have a symbiont load half that of their mother cell. The interplay between horizontal and vertical transmission of symbionts selects a particular symbiont distribution across the host population. We investigate how this distribution responds to different environmental conditions, in particular light levels, and how host control mechanisms may evolve. Note that our model is constructed in reference to the specific relationship between P. bursaria and Chlorella, but is readily reparameterised so as to be applicable to other photosymbioses. Hence we formulate our model using the language of a general symbiotic relationship between a heterotrophic host and a phototrophic symbiont, and parameterise it using data available in the literature on the P. bursaria-Chlorella symbiosis.
The paper is organised as follows. In Section 2 we derive our model. Incorporating nutrient trading and limitation allows us to formulate the host growth rate in terms of symbiont load and nutrient availability, highlighting the different strategies available to the host. This leads to the inclusion of general host control mechanisms, with particular attention paid to their impact on symbiont distribution via horizontal transmission. We then parameterise our model with respect to the specific P. bursaria-Chlorella relationship. In Section 3 we perform numerical simulations of our model, highlighting the roles played by host control and irradiance in determining population equilibria. We then employ adaptive dynamics to illustrate how host control is able to evolve in Section 4. We conclude by discussing our findings and intentions for future investigation in Section 5.

The model
We describe the distribution of symbionts among the host population by defining the time-dependent set of variables φ = (φ 0 , φ 1 , . . .) T , where each entry φ k (t), k ∈ N 0 , of the column vector φ represent the concentration of hosts with k symbionts at time t. We assume that the composition of the population changes due to the following processes: • Cytokinesis of host cells, at rate c k . We assume that the host and symbiont cell cycles are synchronised so as to be concurrent on the appropriate timescale; thus a host containing k symbionts divides into two hosts that each contain k symbionts. Also, we suppose that host cytokinesis is mediated by host population density.
• Death of host cells, at rate d k . Host death is independent of population density.
• Symbiont gain via ingestion of free-living potential symbionts, at rate g k .
• Symbiont loss, at rate l k k, where we assume hosts lose symbionts at a rate proportional to their symbiont load.
Note that the first process encodes vertical transmission of symbionts, while the third and fourth encode horizontal transmission via the free-living population.
Empirical evidence supports the hypothesis of cell cycle synchrony; Takahashi et al. [29] found that symbionts divided only when host cytoplasmic streaming ceased, which occurred just prior to host division. Each symbiont divided approximately once, resulting in two daughter cells with symbiont loads approximately equal to that of the mother before cessation of cytoplasmic streaming. This is opposed to, for example, a timescale-separated situation in which symbiont and host division is not concurrent, or a complete lack of synchrony as expected in an evolutionarily young symbiosis. In addition, we suppose that maintaining symbionts diverts nutrients away from cell growth, so that an excess of symbionts results in a net detrimental effect on the host (parasitism).
We have simplified the gain and loss processes to include only those ingestion events which result in the retention of a new symbiont and assume that loss of symbiont results in ejection from the cell. Although digestion of unwanted symbionts is more likely, for the sake of simplicity we do not explicitly account for this, nor for predation of free-living potential symbionts, in our model. We instead consider such effects to be sufficiently accounted for by the general feeding behaviour of the host (cf. Section 2.1).
Note that all rates depend on the symbiont number k, allowing us to explicitly incorporate the costs and benefits different symbiont loads bring to the host. Moreover, each of the two horizontal transmission processes yields a potential mechanism for host control of the symbiont population, by altering the gain and loss rates in response to the cost/benefit trade-off inherent to the symbiosis.
Taking all this into account, we see that the system varies according to the infinite set of ordinary differential equations where (2) holds for all k ∈ N and we define the total host population K is the carrying capacity of the system. Note that if we define φ −1 ≡ 0, setting k = 0 in (2) yields (1). We write the equation for φ 0 explicitly in order to emphasise the differing behaviour at the boundary k = 0; in particular, that the only way for a host to leave the aposymbiotic state φ 0 is to gain a symbiont via ingestion of an organism from the free-living population.
(1)-(2) can be written using matrix notation as where the entries of the matrix A are given by and those of B by for (k, j) ∈ N 2 0 . Note we have separated (4) into linear and nonlinear parts, with the coefficients of the linear part yielding a tridiagonal matrix, A, and those of the nonlinear part yielding a diagonal matrix, B.
The symbiotic relationship is characterised by a transfer of nutrients from the host to its symbionts in exchange for photosynthetically fixed carbon. Identifying carbon and nitrogen as the primary elements limiting growth, we describe this nutrient trading in the following manner. We assume that the host uses carbon and nitrogen in the ratio λ h , and hence host growth rate is limited by where C k (N k ) is the net intake rate of carbon (nitrogen) used directly by the host for growth, and is dependent upon symbiont number due to nutrient trading. We assume that the uptake of nutrients by phagotrophy (by the host) or photosynthesis (by the symbiont) occurs according to a Holling type II functional response [11,30]. Each symbiont then yields a proportion z s ∈ [0, 1] of its photosynthate to its host; in return, the k symbionts receive a share of the nutrients (in particular, nitrogen) obtained by the host via phagotrophy. We assume that the total nutrition released by the host to its symbionts also varies according to a Holling type II functional response, with maximum z h ∈ [0, 1]; hence, as the symbiont load increases, the host gives up an increasing proportion of its phagotrophically obtained nutrition, but the share received by each individual symbiont decreases. We therefore have where C F and N F are the amount of usable carbon and nitrogen provided by the hosts food, given by and C L is the amount of usable carbon provided by symbiont photosynthesis, given by Here a represents the maximum symbiont photosynthesis rate and b the maximum host phagotrophy rate; F is the amount of bacterial food available to the host and L the light level, each assumed constant; κ F , κ L and κ h are the appropriate half-saturation constants; η L , η C and η N are the conversion efficiencies of photosynthetically obtained carbon, phagotrophically obtained carbon and phagotrophically obtained nitrogen to cell growth; and λ F is the C:N ratio of the hosts food source. Thus the nutrient trading is characterized by the host and symbiont trading traits z h and z s and the half-saturation constant κ h , relating to how a host distributes nutrients among its symbionts. (9)). We rearrange this condition to obtain where compares the λ h -weighted nutrient intake via phagotrophy to the per capita symbiont carbon production (note that photosynthesis yields no nitrogen, and so the denominator contains only a carbon term). From the point of view of the host, µ < 0 represents carbon-deficient food and µ > 0 represents carbon-rich food. Roughly speaking, the closer |µ| is to zero, the more efficient phagotrophy is, in the sense that nutrient intake is closer to the optimal ratio λ h . Equation (13) has at most one non-zero solution k = k λ , where We can discount the negative square root as that solution of (13) is always negative (unless κ h = 0, in which case the negative square root yields k λ = 0; we assume κ h > 0 throughout). As z h , z s ∈ [0, 1], inspection of (15) indicates that k λ is real and positive only if µ < 0, i.e. if food is carbon-deficient. If µ = 0 then k λ = 0, representing the fact that symbionts provide no nutritional benefit when the host feeds on prey with a C:N ratio which is precisely that required for host growth. If µ > 0 or z s = 0 then (13) has no non-negative solutions and (15) is physically meaningless; rather, the best strategy for a host in these cases is to divest itself of symbionts.
From inspection of C k (8) and N k (9), we see that hosts are carbon-limited for 0 ≤ k < k λ , and nitrogen-limited for k > k λ . Hence, if k λ < 0 hosts are nitrogen-limited no matter their symbiont load. Note that k λ is in general not integer-valued, and so a perfect balance of nutrients according to the ratio λ h is in general unobtainable since the symbiont load k is an integer. The best a host can do in practice is choose the symbiont load which minimises C k − λ h N k , i.e. that for which C k − λ h N k is closest to zero, thus minimising any nutrient surplus.
Associated with each symbiont there is also an upkeep cost, separate to nutrition, incorporating such processes as the maintenance of the perialgal vacuole, production of protein transporters, etc. For simplicity, we assume this cost to be proportional to symbiont load. Note that we discount the possibility that symbionts produce any benefits in addition to nutrition, in order to focus our attention on the nutrient trading. In light of the above discussion of nutrient trading and limitation, we therefore write the host cytokinesis rate as where we define the ramp function α c is the rate at which a host converts nutrition into growth (and therefore cytokinesis), while Q is the additional upkeep cost per symbiont. Note that above a certain threshold symbionts become parasites, as the cost of upkeep becomes so burdensome to the host as to outweigh any nutritional benefits. Hence the host death rate is comprising the effects of an excessively burdensome symbiont load and the base death rate in the absence of symbionts, denoted by the constant α d .
We can now define the net host growth rate as Thus we see that a hosts symbiont load directly effects its potential for growth. In light of the analysis leading to the derivation of the optimum symbiont load k λ (15), we have where ⌊·⌋ denotes the integer part of a given real number and k λ is defined in (15). When k λ < 0, hosts are always nitrogen-limited (C k > λ h N k for all k ∈ N 0 ) and r k is a monotonically decreasing function from r 0 . In such a situation symbionts provide no benefit to the host. As our aim in the present article is to model photosymbiosis, we discount this situation as unrepresentative of the biological reality we are interested in, and shall not discuss it further.
Assuming k λ ≥ 0, we can identify three parameter ranges, in each of which the behaviour of r k as a function of k is qualitatively distinct. From (20), we see that r k has a local maximum at k = k max , where and C k and N k are defined in (8)- (9). Furthermore, we can see that the first line of the right-hand side of (20) vanishes at k = k 0 , defined as where k 0 is not necessarily integer-valued. Thus, provided k 0 is real and positive, r k has a local minimum at k = k min , defined as Inspection of (20) indicates that if C L z s < Q then r k is monotonically decreasing; moreover, if C L z s ≥ Q then k min is always real-valued. We shall henceforth assume the latter, as otherwise the per capita cost of maintaining the symbiosis outweighs the gain, in which case the model does not represent a biologically relevant relationship. We can therefore define the three regions in which r k exhibits qualitatively different behaviour as follows: • k min ≥ k max . Here r k+1 > r k for all k ∈ N 0 , yielding a host growth rate which is a monotonically decreasing function of k. Hence symbiont-free hosts exhibit the highest growth rate.
• 0 < k min < k max . In this region k = k min defines a local minimum of r k , with local maxima at k = 0, k max . The existence of two maxima in the growth indicates a potential choice of strategy for the hosts, the precise details of which will depend heavily on horizontal transmission of symbionts.
• k min ≤ 0. r k now has a global maximum at k = k max , with a local minimum at k = 0 and a global minimum at k = ∞. The optimal host strategy is a symbiont load of k max ≈ k λ .
The functional form of this growth rate for the biologically important region k min ≤ 0, and the double-peaked scenario 0 < k min < k max , can both be seen in figure 1.
We now introduce the possibility of a host actively managing its symbionts over and above metabolic provision, as opposed to them being simply a passive load. In effect, the control of its symbiont load becomes a trait of the host. We suppose that control is implemented via the symbiont gain/loss terms in (4), thus enabling hosts to bring about dynamic changes in the symbiont population, and that hosts choose to alter their symbiont load on the basis of increasing fitness, i.e. increasing r k . This mathematical formulation serves as a proxy for the underlying physiological factors influencing the choice between investing in symbionts, or shedding them.
We therefore define and where Θ(x) is the Heaviside step function with Θ(0) = 0 and we assign r −1 = 0. If γ = 0, there is no active host control and horizontal transmission of symbionts is a purely passive process. If γ is non-zero, however, the host manages its symbiont load in order to increase its growth rate, increasing either its gain or loss rate as required in order to achieve this. The terms inside square brackets in (24)-(25) ensure that hosts do not get trapped at a local minimum of r k ; rather, hosts escape the minimum by either gaining or losing a symbiont, with equal proportions choosing each of the two strategies. In contrast, if a host finds itself at a local maximum it reverts to passive symbiont gain and loss only. If its symbiont load should change due to passive processes, host control kicks in once more to restore the optimal state. For simplicity, we assume hosts invest equally in both their gain and their loss rates, increasing each by an equal proportion of the respective passive values. We include the light-dependent factor in g k as a proxy for the free-living population of potential symbionts (cf. the rate of algal photosynthesis per cell in (12)), which we assume constant with respect to time. This ensures that there is no ingestion of new symbionts in the dark, as purely phototrophic organisms cannot survive in the absence of light to drive photosynthesis. Our chosen formulation of the host control via a step increase in an appropriate parameter is probably the simplest from a mathematical point of view; however, it incorporates sufficient biological detail as to yield a useful means by which to investigate the phenomenon of host control. It is clear from the structure of the discrete Burger's equation (4) that horizontal transmission is key in determining the symbiont distribution. This observation motivates the definition of the approximate rate of horizontal Although this is an approximate formula only, it is nonetheless informative as to the qualitative behaviour of the symbiont distribution. Roughly speaking, if v k is positive then symbionts flow to the left, while if v k is negative they flow to the right. Of course, this picture is complicated by interplay between horizontal transmission and host population growth. However, one feature of (26) is especially useful; when v k decreases through zero, horizontal transmission acts to create a net flux towards this point from both directions, in effect creating an attractor for the symbiont dynamics. Thus we expect a peak to form at or near this point, depending on how strongly horizontal transmission is acting.
We plot r k (19) and v k (26) for two qualitatively different parameter regimes in figure 1. In the scenario depicted by the left-hand panels, in which nutrient trading is relatively cheap for the host, r k has a single peak. Thus increasing the host control γ from zero acts to shift k + closer to k max , the point at which r k is maximal, and at the same time increases the magnitude of the rate of horizontal transmission. On the other hand, the right-hand panels depict nutrient trading which is relatively expensive for the host, yielding a growth rate with two maxima. Increasing γ from zero in this instance initially has the same qualitative effect as before. However, once γ increases over a certain threshold then a second point appears at which v k decreases through zero. Thus we now have two possible attractors, separated by a repellor at which v k increases through zero. We discuss the consequences of this further in Section 3.
We denote the possible attractors, i.e. the points at which v k decreases through zero, by k + . These can be found analytically. To this end, we define to be the four possible (and not necessarily integer-valued) solutions of v k = 0. The subscripts indicate which of the symbiont gain (first subscript) and loss (second subscript) processes are merely passive (subscript 0) or under host control (subscript 1) in the region of k-space that particular solution falls. Then k + can be determined by calculating whether v k is positive to the left and negative to the right of each of these solutions. Note we assume β g and β l are strictly positive throughout. We note first that if γ = 0 then k + = k 00 . For γ > 0, the situation is more complicated due to our formulation of the gain and loss rates (24)- (25) in terms of the discontinuous Heaviside function. Assuming a single-peaked growth rate, so that k 0 < 0 (cf. (22)), we have see the left-hand panels of figure 1. If instead k 0 ≥ 0, yielding a growth rate with two maxima, then see the right-hand panels of 1. Note that k min is never an attractor if γ > 0, and neither is zero if β g > 0; hence we have omitted these possibilities from (29). We have k 01 ≤ k 00 ≤ k 10 , for all physically realistic parameter values. We can therefore interpret the effects of increasing host control as follows. When γ = 0, both (28) and (29) simplify to k + = max(k 00 , 0), i.e. horizontal transfer is strictly passive. As γ increases from zero, k 10 increases while k 01 decreases. Thus k + approaches the local maxima of r k as γ increases, with k + double-valued if 0 < k min < k max . We note there is a third possibility, in which k = 0 is a global maximum; in this case increasing host control pushes k + closer to zero, in effect destroying the symbiosis. We dismiss this case as being biologically uninteresting, and mention it only for completeness. We can therefore see that there are four possible outcomes of the population dynamics: • Stable symbiosis. k + is single-valued; the peak drifts until it reaches k + , where it remains.
• Dichotomous symbiosis. k + is double-valued. Hosts may potentially split between the two strategies of low and high symbiont loads.
• Near loss of symbionts. This corresponds to the limit k + → 0 as β g /β l → 0. In practice, a large but finite β l results in a very low but non-vanishing symbiont load.
• Extinction due to symbiont overload. This corresponds to the limit k + → 0 as β g /β l → ∞. Gain of symbionts is sufficiently higher than loss to push the population into the region of negative growth (r k < 0), at which point the cost of maintaining symbionts is too much and the hosts die, leading to host population crash.
We shall focus on the first two outcomes as biologically relevant, with particular attention paid to the stable, single-peaked symbiosis.
We have estimated the model parameters from data available in the literature. We summarise the values used in table 1, and outline the rationale behind each choice below. Some data have been converted into the units used in the present work; we omit the details where this is trivial, e.g. hours to days. All parameters are given to two significant figures.
• C:N ratios. Finlay and Uhlig ( [31], table 1) calculated the elemental composition of various protozoa, of which the closest relative to P. bursaria is P. caudatum, which has a C:N ratio of 3.5 by weight, or 4.1 by amount of substance. They also found the C:N ratio by weight of a freshwater bacterial sample to be 3.9, or 4.6 by amount of substance. This is in agreement with the value given by Fagerbakke et al ( [32], table 2) for various native aquatic and cultured bacteria.
• Conversion efficiencies. Herrig and Falkowski [33] found the conversion efficiency of photosynthetically produced carbon into growth by algae to be 85%, which we shall use for η L . Jones et al [34] found that nitrogen growth efficiency was around 2.5 times the carbon growth efficiency of a marine copepod when nitrogen was in short supply, and so intake is maximised at the expense of efficiency. As hosts must uptake enough nitrogen for themselves and their symbionts, we expect nitrogen intake to be similarly prioritised in our system. Furthermore, Mauclaire et al. [35] found that 12% of predated bacterial carbon was transformed into protist biomass in a bacteria-protist predator-prey system. These data lead us to assign the values η C = 0.12 and η N = 0.3 to the phagotrophic conversion efficiencies.
• Photosynthesis. In [12], Finlay et al. analysed symbiont photosynthesis in a naturally occurring population of mixotrophic ciliates. In particular, they fitted photosynthetic data to a Holling type II function as in (12) • Phagotrophy. Fenchel [11] performed a comprehensive investigation of suspension feeding in ciliates, including data from the two Paramecium species P. caudatum and P. trichium, which we shall use in lieu of data on P. bursaria. Figure 2 of [11] provides an estimate for the maximum ingestion rate as 4.32 × 10 5 µm 3 d −1 . Table 2 in [32] yields the estimated values for the cellular C:volume ratio of 8.3 × 10 −9 µmol µm −3 . As the C:N ratio is λ F = 14/3, the proportion of food by amount of substance which is carbon is Figure 6 in [11] provides a value of 4 × 10 7 µm 3 ml − 1 for the half saturation constant. We estimate the volume of a bacterial cell to be 2µm 3 from the data in table 1 of [32], yielding κ F = 2 × 10 7 cell ml −1 . We take the bacterial food concentration to be 5.8 × 10 6 cell ml −1 .
• Host growth. The precise values of α c , α d and Q are difficult to determine, as we have formulated our model in terms of a growth rate r k dependent on individual symbiont loads, whereas data is calculated at the population level. However, the data of Karakashian [9], who carried out a series of growth experiments in differing environmental conditions, yields enough information for us to make reasonable estimates. Bleached P. bursaria grown in sterile media were found to simply die due to the complete lack of available nutrition, providing an estimate for the base death rate α d = 0.7 d −1 ([9], figure 1). Green P. bursaria grown with abundant food and light had a mean daily fission rate of around 1.9 d −1 ([9], table 2). Furthermore, we know from the analysis around (22) that for symbiosis to be desirable the upkeep cost Q per symbiont must be less than the nutritional benefit per symbiont. Combining these considerations leads us to choose α c = 1.8 × 10 4 cell µmol nutrient −1 and Q = 5 × 10 −8 µmol nutrient −1 cell −1 d −1 . We choose K = 1000 cell ml −1 as representative of a typical carrying capacity.
• Nutrient trading. z h , z s and κ h are traits we vary in order to investigate the symbiosis, and so we shall not fix them at any particular value.
• Horizontal transmission. In a reinfection experiment [13], Karakashian reported that aposymbiotic P. bursaria exposed to a highly concentrated algal population ingested a mean of 285 in 2.5 mins, and after 23.5 hours a mean of 60 of these remained. We shall therefore take 60d −1 as the maximum ingestion rate, as we assume that if the hosts had means and motive to ingest more symbionts then more than 60 would have remained. We shall assume that, as experimental irradiance was high, host control was in effect, and so the passive ingestion rate β g is lower than 60 d −1 ; we choose β g = 30 d −1 and vary γ from zero. Data appears to be scarce to nonexistent for rates of symbiont loss. We therefore simply assume a low passive loss of β l = 0.1d −1 .

Population equilibria
We now present numerical solutions of the model (4), highlighting the key ecological processes leading to population equilibria. We truncate the system at k = k trunc , where k trunc is chosen so that φ ktrunc is exponentially small and therefore has a negligible effect on the solution near the symbiont peak. For the present purposes, setting k trunc = 500 turns out to be sufficient; although the growth rate r 500 is usually positive for the parameter values in table 1, horizontal transmission is sufficiently negative as to render population equilibria in which φ 500 ≈ 0. We calculate steady-state solutions of the nonlinear problem (4) using Newton-Raphson iteration. We also calculate In figure 2 we present steady state solutions of (4) with a single-peaked growth rate distribution, corresponding to the left-hand panels of figure 1. We see that the solution is closely approximated by the dominant eigenvalue of the associated linear problem. In fact, our numerical investigations showed that this agreement is excellent for a wide range of relevant parameter values, suggesting horizontal transmission rapidly organises the symbiont distribution. Furthermore, the location of the peak is approximately given by k + (not shown for the sake of clarity; cf. (28)), indicating that v k (26) is an excellent description of the rate of horizontal transmission. We can also see how increasing host control strength γ moves the symbiont peak closer to the optimal value k λ , and reduces its variance, indicating that hosts are maximising their benefit from the symbiotic relationship as γ increases.
In figure 3 we present steady-state solutions for a double-peaked growth rate distribution, corresponding to the right-hand panels of figure 1, although we note that the per capita cost Q in figure 3 is double that of figure 1. Again, agreement between steady state solutions and dominant eigenvectors is excellent, although we omit the eigenvectors for clarity. For low γ, the peak is near k λ ; as γ increases, the solution becomes double-peaked before losing the upper peak altogether as γ increases further. In this case, increasing γ without bound will have the effect of destroying the symbiosis, as increased host control moves the lower peak ever closer to zero. We note that not all double-peaked growth rates, or even double-valued k + , yield corresponding double-peaked symbiont distributions. The deciding factor appears to be the net flux towards an attractor; the upper peak is usually at an advantage in this respect as hosts originating to the right of k max are attracted to the upper peak and hence cannot reach the lower.
A key environmental property affecting the symbiosis is light. In figure 4 we plot the total host population Φ and the mean symbiont load of steady state solutions to (4), for four levels of host control. We note that in each case the agreement between the numerically computed values ofk and the approximate prediction k + is excellent, but neglect to plot this information to preserve clarity. We see in figure 4b that when there is no host control, the mean symbiont number simply increases with light according to the functional response incorporated into g k (cf. (24)). With a moderate level of host control, the mean symbiont number increases from zero to a local maximum, at which point it decreases in line with k λ before increasing again as passive gain overcomes host control due to an increasing free-living population. Strong enough host control prevents this, ensuring optimal nutrient intake at all light levels above a threshold, correlating with the observed increase in host population with higher levels of host control. We note that there is an initial dip in host population as irradiance increases from zero, indicating that at first symbionts are not photosynthesising enough to provide a net benefit to their hosts. Increased host control reduces this dip as the hosts invest in expelling unwanted symbionts. Note the sharp increase ink at low light; the solution changes qualitatively at this point, with a discontinuous change in the location of the peak.

Adaptive dynamics
Having discussed the ecological aspects of the model in detail, we now turn our attention to the evolution of the symbiosis. In particular, we shall investigate the evolution of host control by allowing the strength of control, γ, to evolve in an adaptive dynamics framework. The benefits of increased host control are likely to be balanced by an associated cost, which we model by taking the per capita cost of symbiosis maintenance to be a monotonically increasing function of γ, specifically where the Q i , i = 1, 2, 3, are constants. Thus the cost of symbiosis increases from an initial value of Q 0 at γ = 0 to a maximum of Q 0 + Q 1 at γ = ∞, with Q 2 determining the rate at which this transition occurs. This limiting behaviour of the per symbiont cost provides a trade-off curvature as the trait evolves. The adaptive dynamics approach assumes a separation of timescales between ecological and evolutionary dynamics, in that the ecological equilibrium is continually updated in evolutionary time by a series of successful invasions by initially rare mutants. A host with a mutant trait γ ′ must be able to overcome the background state provided by the resident steady state population with trait γ. Hence the initial dynamics of a newly-appeared mutant ψ(t) are governed by where Φ * = ∞ j=0 φ * j and the φ * j are a steady-state solution of (4) in the absence of mutants, i.e.
Note that unprimed rate functions refer to the resident population with trait γ, and primed variables refer to the mutant population with trait γ ′ . We can see that (33) is an eigenvalue problem. Therefore, if there exists an eigenvalue with positive real part corresponding to an eigenvector with no negative elements, the initially rare mutant population will exhibit exponential growth until nonlinear effects become significant. We shall make the usual adaptive dynamics assumption that such an invasion is always successful, with mutants replacing residents and achieving a new equilibrium, distinct from φ * due to the updated trait value γ ′ .
The eigenvalues of (33) must be calculated numerically. To this end, we first found the solution curve comprising solutions to (34) for different values of the resident trait γ across the desired range. Then, for each resident trait we calculated the eigenvalues of (33), for each value of the mutant trait γ ′ across the same range. Following the adaptive dynamics methodology, we differentiate between the regions in which at least one eigenvalue corresponding to an eigenvector with no negative components has positive real part, and those in which none do. The former regions of trait space are therefore those in which a successful invasion may occur. We plot our results in figure 5 for various values of Q 1 , taking Q 0 = 2 × 10 −8 and Q 2 = 0.5. We see that if the cost increases relatively gradually with γ, a moderate level of host control will evolve. As Q 1 increases, an evolutionary steady state which is not convergence stable emerges from the origin, resulting in a barrier to the evolution of host control; if a mutation appears which is large enough, however, there still exists an evolutionary stable strategy with γ = 0. As Q 1 increase further, the two evolutionary steady states approach one another and coalesce before vanishing, in which case evolution favours no host control. We note that evolutionary stable strategies are associated with host population maxima, and evolutionary unstable strategies with population minima.

Discussion
We have constructed a deterministic model of the endosymbiosis between a heterotrophic host and a phototrophic symbiont. By formulating the host growth rate in terms of its nutritional state, we were able to explicitly incorporate nutrient trading between a host and its symbionts into our model, allowing us to investigate the consequences of symbiosis for host fitness. A novel feature of our model is an explicit mechanism by which the host can exert a degree of control over the number of endosymbionts via adjustment of the rates of loss and gain of symbionts. The strength of this host control is captured by a single parameter in our model.
Host control is an essential addition to a model of this symbiosis in order to describe the behaviour of holobiont in response to variable levels of light. The model presented here predicts that the optimal symbiont load k λ (15) decreases monotonically with irradiance from infinity at zero light. This is because fewer symbionts are needed for the same gain as photosynthetic output increases, but this gain can only be realised if the host divests itself of excess symbionts, thus diverting nitrogen from symbiont nutrition back to its own growth. Such a response requires a level of host control over its symbiont load-without it, the symbiont population will increase with light to the detriment of the host. Thus, as light increases from zero, the symbiotic relationship moves from parasitism through commensalism to exploitation of one party by the other, depending on the level of control exerted by the host.
The model leads to the following predictions at different light levels. When environmental irradiance increases from zero, after an initial dip in which symbiosis is slightly detrimental to the host, we observe a sharp increase in symbiont load. This is mainly due to an increase in ingestion as the free-living population of algae increases from zero in the dark. As light levels continue to increase, at a certain irradiance the mean symbiont load coincides with the optimal load k λ . According to the level of explicit control imposed by the host on the symbiont, the model predicts that either the symbiont load will continue to increase, proportional to the rate of photosynthesis, or will remain at k λ and therefore decrease with increasing light. If host control is strong enough, this decrease is monotonic and the mean symbiont load remains approximately equal to k λ . However, for weaker levels of host control, the symbiont load once again begins to increase due to the free-living population becoming too large for the host to effectively manage (cf. figure 4). Thus the differing strategies at higher light levels introduce conflict between host and symbionts, with the possibility that symbionts may counter-adapt defences against host control. Our model therefore demonstrates that adaptation to different light levels may not lead to easily predictable results, as coevolution may drive the symbiosis to differing regimes. In particular, strains adapted to higher light levels should exhibit a greater level of host control in order to overcome the increased symbiont gain due to increased growth of free-living algae and prevent exploitation of the hosts by their symbionts.
This pattern of behaviour is precisely seen in recent experiments on this host-symbiont relationship. Lowe et al. [36] show that a peak in symbiont load is observed at a relatively low irradiance, and the symbiont load then decreases monotonically as the light levels increase. In our model this suggests that the symbiotic relationship of the strain studied has evolved to have a high level of host control. To confirm that such an evolutionary end point is a likely result of our model we have presented an adaptive dynamics approach using the host control strength γ as the evolvable trait. The control is traded off against the per symbiont cost to the host, representing the investment of the host in explicit mechanisms to restrain the symbiontnutrient consumption is explicitly accounted for elsewhere in the model. Our evolutionary model demonstrates that when the cost to the host of imposing control is low then it is favourable to adopt a high level of host control; there is a convergent and evolutionary stable state for a finite value of γ. However, the fine detail of this evolutionary end point is sensitive to the per-capita costs of control, which suggests that different strains are likely to have evolved different levels of host control and this should be reflected in their response to a light gradient.
If the per capita cost Q increases too severely with the evolution of control then the convergent and evolutionary stable state is γ = 0, i.e. hosts do not evolve control. In this parameter regime our adaptive dynamics description no longer encapsulates the essential evolutionary pressures, with evolution beginning to act on the passive gain and loss rates via the parameters β g (24) and β l (25), for example. A likely end result of such evolution is the transition back into a regime in which conditions are favourable for the evolution of host control.
There exists also a more drastic scenario which we have not covered in Section 4. If Q increases sufficiently it can cause a transition from the solution branch in which the symbiont peak is near k λ to that with a peak to the left of k 0 . In this case, increasing host control acts to destroy the symbiosis as the symbiont peak moves towards k = 0 (cf. the discussion around figure 3 in Section 3). Furthermore, γ is then free to increase without bound until the symbiosis is completely destroyed. However, once again we find our adaptive dynamics description failing to capture all relevant evolutionary effects; in this case, we expect ever increasing investment in host control to begin to have a detrimental effect on host growth in other ways than only increasing Q, for example decreasing the cytokinesis rate α c (16) or increasing the death rate α d (18).
Double-peaked growth rate distributions, in which the host has a choice of strategies, occur when nutrient trading is expensive for the host, thus reducing its benefit from the symbiosis. In reinfection experiments [13,37], it is often observed that the resulting symbiont load distribution is doublepeaked. Our model suggests that this could be due to a decreased release of photosynthate from new symbionts to their hosts, potentially a result of the time spent living autonomously, or perhaps due to a delay between ingestion of symbionts and commencement of nutrient trading. Moreover, our results indicate it is likely that one peak is transient, and given enough time, will decay leaving only a single peak. As discussed above, the double-peaked growth rate distribution provides a mechanism, as Q increases, whereby it is possible for the evolutionary system to select a non-symbiotic state, but without explicit modelling of the symbiont dynamics and ecology this remains a speculation.
The construction of a more realistic mechanism by which cell-cycles become synchronised, by incorporating nutritional dependence into symbiotic growth rates, is a desirable feature currently lacking in our model. This also permits the inclusion of other host control mechanisms, implicit here, whereby the host manages the nutrient supply to its symbionts in order to maintain a stable population, perhaps including such dynamic nutrient trading as that described in [28]. The precise nature of cell-cycle synchronisation is at present poorly understood, and including such effects in our more comprehensive model would build greatly on the simpler approaches attempted here and previously [20,21,22,28]. We note that, in contrast to [28], for example, our model highlights the importance of horizontal transmission in determining the stable symbiont distribution. In fact horizontal transmission of symbionts is necessary, when combined with host control, in order to react appropriately to environmental changes (cf. [36]).
The model presented here forms an excellent basis for further study of endosymbiotic relationships. Our model gives a number of detailed predictions regarding the photosymbiotic relationship between P. bursaria and Chlorella, in particular, the dependence of symbiont load on light and host control. These predictions provide insight into the underlying mechanisms and trading relationships between the two partners. Subsequent investigations into these details may permit greater understanding of the nature of the relationship between the two partners including its equality, its history, its future trajectory and ultimately greater understanding of this important route to complex lifeforms.