Lotka-Volterra, replicator dynamics, and stag hunting bacteria

Last time in the Petri dish, I considered the replicator dynamics between type-A and type-B cells abstractly. In the comments, Arne Traulsen pointed me to Li et al. (2015):

We have attempted something similar in spirit with bacteria. Looking at frequencies alone, it looked like coordination. But taking into account growth led to different conclusions […] In that case, things were more subtle than anticipated…

So following their spirit, I will get more concrete in this post and replace type-A by Curvibacter sp. AEP13 and type-B by Duganella sp. C1.2 — two bacteria that help fresh water Hydra avoid fungal infection. And I will also show how to extend our replicator dynamics with growth and changing cell density.

Although I try to follow Arne’s work very closely, I had not read Li et al. (2015) before, so I scheduled it for a reading group this past Friday. I really enjoyed the experiments that they conducted, but I don’t agree with their interpretations that taking growth into account leads to a different conclusion. In this post, I will sketch how they measured their experimental system and then provide a replicator equation representation of the Lotka-Volterra model they use to interpret their results. From this, we’ll be able to conclude that C and D are playing the Stag Hunt — or coordination, or assurance, pick your favorite terminology — game.

Measurements in the experimental model

The challenge of working with microscopic organisms is that there are a lot of them, you can’t easily see them, and so they are hard to count. The two Betaproteobacteria Curvibacter sp. AEP1.3 and Duganella sp. C1.2, that Li et al. (2015) selected, offer a way around this. When placed in very low concentrations on R2A agar plates, the cells grow into colonies that become visible to the naked eye. More importantly, the two different cell-types grow into colonies at different rates, with D colonies becoming visible to the naked eye after 2 days, and C colonies taking 4 days. Due to the heavy dilution, each colony could be assumed to have started with a single cell, and so a count of the colonies provided a count of the cells. For each measurement, Li et al. (2015) did two dilutions, one with a higher factor than the other. If figure 1 is representative then it looks like the colony counts were on the order for ~500 cells for dilution 1 and ~20 cells for dilution 2.[1]

The dynamics of the experimental system were run in a test tube with the cells in solution. This has the benefit of working uncontroversially with the typical inviscid story of replicator dynamics, but carries the cost — for my particular interests of applications in mathematical oncology — of not generalizing easily to cancer cells, which prefer to grow while adhering to a surface. The test tube of growth medium was seeded with an initial concentration of colony-forming units per milliliter, with 14 different initial proportions of C versus D cells. Three times a day for three days, Cleo Pietschke (1) measured the optical density — at a wavelength of 600 nm (OD600) — of cells in the test tube, (2) picked dillution factors based on that measurement, (3) took a small sample of the cells in solution, and (4) plated them for measurement as described in the previous paragraph, with the dilution factor from (2).

Replicator dynamics for growing populations

Let me resume where I left off last time with the replicator dynamics for the proportion p of C interacting with D cells given by:[2]

The easiest way to tell a micro-dynamic story for the replicator equation is with fixed population sizes.[3] But proportions aren’t the only information we have from an experiment, and population sizes aren’t fixed. In fact, the biggest difference between ecological modeling in micro- vs. macro- organisms is that macro-organisms seldom have the opportunity to undergo exponential growth; they are almost always at carrying capacity. So what if we want to model this ecological difference — the fact that the total population grows or that the cell density in the Petri dish changes?

The easiest way to introduce population density into the replicator dynamics, is just to suppose that the whole population follows logistic growth at a rate equal to the average fitness ; or in equation form:

Note that this equation is coupled to the dynamics for p, but the dynamics for p are independent of unless or depend explicitly on . For me, this means that if we want to understand the dynamics of the cell proportions, we gain nothing from knowing the dynamics of the cell density. Explicitly modeling cell density allows us to understand cell density, but it does not change the dynamics of cell proportions. It is still the same game.[4]

But the above equations are equivalent to what Li et al. (2015) consider. We can see this by rewriting the total density in terms of the density of C () and D () types:[5]

and follow the same procedure to get a similar equation for density of type-D cells:

Now, sometimes we might care about the actual counts of the cells, instead of their densities. In that case, suppose that our populations have a carrying capacity of K then we get:

And this is the same as equations 3.1 in Li et al. (2015) in the case that the rest of the paper considers.[6]

Being able to switch between these two equivalent representations makes analysis much easier.[7] For example, to show that if then the dynamics lead towards a ‘stable manifold’ with , Li et al. (2015) have to explicitly consider the time derivative of the distance function in equation 3.4. However in the replicator dynamic representation, you can just read that right off the equation for — remembering that . In fact, you could go further to classify all possible population density dynamics for any — even if they are not strictly positive — in terms of the roots and sign-changes of .

Finally, since this representation shows that the proportion dynamics are not dependent on the density, it means that we can just look at the proportion dynamics in figure 2 — reproduced above — of Li et al. (2015). Since these dynamics bifurcate towards all-C or all-D, depending on if they start with a proportion of C higher or lower than p*,[8], I would conclude — as the authors suggest — that this is a stag hunt — or coordination, or assurance — game.

Quadratic fitness function for Duganella sp. C1.2

Of course, looking at the qualitative features of the dynamic profile is not the only — and not the best — way to figure out game dynamics. A more quantitative approach would be to actually measure the fitness function.[9]. Since, Li et al. (2015) kept their bacteria in the growth-phase during experiment, they estimated fitness by plotting the log of each of the two cell type’s population and then taking the slope of the line of best fit as their estimate of the growth rate.[10] In figure 5b — reproduced below — they plotted these estimated growth rates versus initial frequency of C.

From this figure, they then used both Akaike and Bayesian information criteria to select the curves of best fit from among linear, quadratic, and cubic models as the empirical fitness functions:

These are plotted above as dashed lines in the figure above. As you can see from the figure and the fit, is not linear in p, but quadratic. This is the most exciting conclusion of the work for me. It reinforces my belief that evolutionary game theorists should expand their heuristic models from linear to higher order gain functions, even if we don’t have good micro-dynamical stories for this transition.[11] With this work, Li et al. (2015) provide a great contribution to the nascent literature of closely-coupled empirical measurements and mathematical models of microscopic dynamics in evolutionary game theory. A move that I hope to expand on in future work and in discussions here on TheEGG.

Notes and References

To me, this seems like a very noisy measurement and I would like to have at least seen error bars on figures 2 and 5, if not a full discussion of metrology.

If you are following along with Li et al. (2015) in hand then note that I use p and w where they use x and r.

Unfortunately, this easiest story seems to have led many to conclude that replicator dynamics is only for constant (or nearly constant) populations, and to decry it as inappropriate for growing populations. I don’t understand this sentiment.

To make density dependence interesting for me, there would have to be a case where something different happens to the density of the cells depending on the initial densities or proportions of cells. As is the case in the ecological public goods game (Hauert et al., 2006; 2008).

Notice that although the first equation is interpreted as replicator dynamics and the second is logistic growth, they both have the same form. If we wanted then we could think of them abstractly as a factoring of 3-strategy replicator dynamics between density of A-types, density of B-types, and density of ‘free space’-types, where .

The equations do differ in the case, with my approach giving us:

while Li et al. (2015) have the a instead of the in the first line, and vice versa in the second. Li et al. (2015) present their equations 3.1 as just the standard competitive Lotka-Volterra equation, but the competitive Lotka-Volterra equation allows for arbitrary couplings and . Their equation would be equivalent to setting and , while my version would have and . Their version seems simpler, but I don’t know of a grounded justification for it or for my version. It certainly doesn’t seem unreasonable — as in my version — to have the coupling depend on both the fitness and the carrying capacity.

Further, if we take units seriously and do our stoichiometry then I think that strengthens the case for my version of equations 3.1. If we let C be measured in [cf units of C] (and D in [cfu of D]), in units of [cfu of C]/[unit of time] (and accordingly for D), and the growth rates as 1/[unit time] then the units of my above equation check-out. But with this choice of units — which is by no means unique, and maybe the authors had a different parametrization in mind — the equations 3.1 from Li et al. (2015) do not type-check:

where I bolded the discrepancy. To overcome this discrepancy, we need to either say that [cfu of C] and [cfu of D] are the same thing or add a conversion factor between C and D to the equations. We could say that this conversion factor is 1, but it seems more natural to recognize that C and D have different carrying capacities and and thus take up differing amounts of resource-space. If we pursue the latter route then the natural conversion factor from D to C is , and appending that factor to the summands in question from equations 3.1 would then yield my equations from the start of this footnote.

If I am mistaken in my reasoning here, or if there is a better choice of units, then I hope that you guys can help explain this to me in the comments.

Note that my point here is not just that 2-population Lotka-Volterra equations can be transformed into 3-strategy replicator dynamics. This trick is well known, and by itself does not establish my main point — that this game can be studied just as a 2-strategy replicator dynamics between C and D — because usually you would end up with a 3-strategy game where (1) the dynamics between your two ‘real’ strategies depend on the dynamics of the fictitious ‘free-space’ strategy, and (2) the proportions of the two ‘real’ strategies won’t exactly correspond to the proportions of their corresponding Lotka-Volterra species. With the transformation in this post, the proportion p corresponds to the proportion of C in the population — or in the case, the proportion of C counted by carrying-capacity use; i.e. — and the dynamics for p are independent from the dynamics for . This might be due to the nature of the transformation that I use. The classic literature (Hofbauer, 1981; Bomze, 1983; Hofbauer & Sigmund, 1998) uses the transformation , while in this post I use for the (i.e. density) representation, or equivalently for the factored representation.

Although it would have been nice to see a more detailed analysis of the crucial around the unstable fixed point. Similar for the stable fixed point that Li et al. (2015) expect from their quadratic fitness fit, since from figure 2 it is hard for me to tell if the population is moving towards all-D, or a stable fixed point very near all-D. Figure 6 gives more evidence for this, but reveals another weakness: the population was not grown for a very long time under K-selection and so we cannot see that the unstable fixed point near all-C in C-D space is in fact unstable, and same for the stability of the fixed point near all-D in C-D space.

Previously, I described how to measure the gain function by looking at the change in the log-odds of strategy C for varying initial proportions of strategy C. Unfortunately, the data in figure 2 is too noisy for this approach; especially without a quantification of the errors of individual measurements.

I might be missing some subtlety of their procedure for measuring the growth rates they use as the y-axis in figure 5b. However, if it is as I describe it in the main text then there seems to be an inherent tension. On the one hand, the authors expect the fitness to change as a function of the proportion of C cells, that is why they plot it versus initial frequency. On the other hand, they measure a single fitness (or growth rate) by fitting a line to the whole 3-day growth curve; but the proportions of cells change drastically over those 3 days (as evident from figure 2). So in the first case, we have to assume that proportions matter, but in the second case, we have to assume that a single growth rate is enough to describe them.

On a final minor point about fitting the data. As Li et al. (2015) describe at the start of section 2.2, they calculated the goodness of fit for these growth rate estimates, which can be a great measure of the error on their estimates of growth rate. I would be interested to see figure 5b with these goodness-of-fit as error bars.

This isn’t to suggest that there are no good micro-dynamical accounts for nonlinear gain-functions. One can look to multiplayer games (Gokhale & Traulsen, 2014; Ohtsuki, 2014) — as the authors discuss — or their close relative the non-linear public goods game (Archetti, 2013; 2014) for inspiration, where the degree of the gain-function is bounded by the number of interaction partners — or one less than the number of players. However, what makes these recent papers exceptional is that they are not the standard mindset in EGT, but one that we have to spread.

Related

About Artem KaznatcheevFrom the Department of Computer Science at Oxford University and Department of Translational Hematology & Oncology Research at Cleveland Clinic, I marvel at the world through algorithmic lenses. My mind is drawn to evolutionary dynamics, theoretical computer science, mathematical oncology, computational learning theory, and philosophy of science. Previously I was at the Department of Integrated Mathematical Oncology at Moffitt Cancer Center, and the School of Computer Science and Department of Psychology at McGill University. In a past life, I worried about quantum queries at the Institute for Quantum Computing and Department of Combinatorics & Optimization at University of Waterloo and as a visitor to the Centre for Quantum Technologies at National University of Singapore. Meander with me on Google+ and Twitter.

If you think of payoffs as arising from interactions (either physical contact or the production of some diffusible factor) then they should depend not only on the fractions of the strategies, but also on the density. Think of an infinitesimal population density: the cells/players will in this limit not interact at all, and hence they should reproduce with some basal rate independent of interactions. As the density grows interactions will become more and more common and contribute more to the reproductive rate.

You are correct, fitness functions could definitely depend on the densities under certain microdynamical assumptions; I don’t think anybody is disputing that possibility. If you choose to make a different set of assumptions then it might depend on the proportion, or even on the ratio of the two populations (for example, see Arditi & Ginzburg, 1989; Akcakaya et al., 1995). The point of my recent posts on operationalization is to get away from armchair theories based on our physical intuitions — especially to get away from my own physical intuitions that I don’t think I have any reason to believe to be in line with biology — and think about how we will measure these functions in the lab.

In the case of checking for density-dependence, there is a clear control to test: seed your experiments with differing densities but same proportions. Li et al. (2015) do this for monoculture — where they show fitness to be independent of density — but don’t report different density seeds for coculture. I don’t know if they performed these controls or not for coculture, maybe one of the authors can chime in and clarify this. In experiments that we plan to carry out here, I will definitely advocate for doing these controls. However, since the point of my post was to engage with Li et al. (2015) on the terms of their paper, I considered — as they do — the case where are functions of the proportion p. If are also functions of then — as I mention in passing in the post — the dynamics of p will depend on the dynamics of , and one will have to explicitly consider the fictitious free-space third strategy.

As for what exact microdynamical story yields explicit dependence — and there definitely are some that will — I am not sure if yours is sufficient. In particular, the Lotka-Volterra representation of the equations (i.e. dynamics for C and D) contain quadratic terms of and so can account — at least to some extent — for increased frequency of interactions with increased density.

If you want to expand more on your point than this comment thread affords then I am always eager for new guest posts from you.

Thanks for your detailed answer. If you are trying to avoid armchair theories then why do you start with the replicator equation? Frequency-dependence is a rather strong assumption, and to me it makes more sense to start with a general model of the type:

dC/dt = f(C,D,t)
dD/dt = g(C,D,t)

and then try to infer the function f and g from experimental data. It might be that f and/or g can be understood best in terms of frequency and/or density, but without any knowledge about the microdynamics it seems best to avoid any a priori assumptions. Frequency-dependence might seem probable from an EGT perspective, but empirical evidence in oncology is basically non-existent.