Evolutionary game theory for physical and biological scientists. I. Training and validating population dynamics equations

Failure to understand evolutionary dynamics has been hypothesized as limiting our ability to control biological systems. An increasing awareness of similarities between macroscopic ecosystems and cellular tissues has inspired optimism that game theory will provide insights into the progression and control of cancer. To realize this potential, the ability to compare game theoretic models and experimental measurements of population dynamics should be broadly disseminated. In this tutorial, we present an analysis method that can be used to train parameters in game theoretic dynamics equations, used to validate the resulting equations, and used to make predictions to challenge these equations and to design treatment strategies. The data analysis techniques in this tutorial are adapted from the analysis of reaction kinetics using the method of initial rates taught in undergraduate general chemistry courses. Reliance on computer programming is avoided to encourage the adoption of these methods as routine bench activities.


Curve-fitting for ternary and higher-order ecologies
In this section, we show how the techniques demonstrated using the two-population system in section 3 of the main manuscript can be extended to a co-culture of three or more subpopulations. We slightly generalize equations 3.1 and 3.2 from the main manuscript by allowing the co-culture to contain an arbitrary number of subpopulations labeled x, y, . . . , z, where we take a notational shortcut by allowing letters of the alphabet, perhaps, to exist between y and z (nonetheless, to simplify figures, most of the discussion in this section uses the example of a 3-population system). We also rename the coefficients A, B, C, and D using indexed coefficients A xx , A xy , etc. Again, it is customary to refer to the rate coefficients as fitnesses. For example, the fitness of population z is labeled f z and equals A zx p x + A zy p y + . . . + A zz p z . Figure 5 illustrates that the parameter values in equations S.11-S.13 are trained using the same method that was used to train the parameter values in equations 3.1 and 3.2. Panel (a) shows that co-cultures of three cell types (x, y, and z) can be prepared so as to be almost purely enriched in cells of type x, almost purely enriched in cells of type y, or, as another alternative, almost purely enriched in cells of type z. The highlighted example is initially dominated by cells of type x. Setting p x ~ 1 and the remaining population fractions p y and p z approximately equal to zero, equation S.12 becomes

Training
In the same way that was done using equation 3.6 in the main manuscript, the coefficient A yx is estimated by measuring the slope of the line tangent to the plot of population y as a function of time t.

S.15
The remaining parameters A xx . . . A zz can be obtained in a similar way, and their values for this example can be found in the following section of this supplemental document.
The final step in training equations S.11-S.13 is to draw a velocity field. Unfortunately, when describing the dynamics of more than two populations, it is difficult to draw arrows using coordinate axes directly analogous to the coordinate axes in Figure 3(h) from the main manuscript. Figure 5 with analogous relationships for the remaining population fractions p y . . . p z . Substituting the trained values of the parameters A xx . . . A zz and the population fractions p x = 0.4, p y = 0.5, and p z = 0.1 into equation S.16 allows us to estimate the change in population fraction p x in the same way that we estimated the change in absolute population x in equations 3.9 and 3.10 in the main manuscript.
Calculating changes in p y and p z then allows us to specify the population composition at the tip of the highlighted arrow in panel (b). In this diagram, each arrow represents a change in population composition corresponding to a time interval of 2 days.

Validation
Once we have sketched a velocity field on the simplex, we can compare it with the trajectory traced out by a dataset that was not initially used to train the parameters in equations S.11-S.13. In the simulation in Before we consider panel (d) a successful validation, it is important to note that a simplex only displays population fractions, not absolute population sizes. This means that the dynamics of the absolute population sizes in (c) must be spot-checked to ensure that the rates of change of absolute population sizes are consistent with equations S.11-S.13. As an example, an approximate tangent line is drawn through the point y = 239 cells at t = 24 days. The slope of this dashed line is approximately +80 cells/6 days, or 13 cells/day. Substituting y = 239 cells, the population fractions p x = 0.33, p y = 0.59, and p z = 0.08, along with the parameters A xx . . . A zz , into equation S.12 predicts that population y achieves a rate of change of dy/dt = 15 cells/day, similar to the 13 cells/day measured. Analogous comparisons would need to be made at a variety of time points in (c) for multiple populations for us to be confident of the empirical accuracy of equations S.11-S.13. These spot-checks can be tedious. Unfortunately, these calculations become especially necessary when studying co-cultures of more than 3 populations. In such cases, validation is not convenient to perform on a simplex, which becomes a multidimensional structure that is difficult to draw on a sheet of paper.

Prediction
Now that we have performed a preliminary validation by using the velocity field in (d) and by measuring slope, as in (c), we can plan additional experiments to further test equations S.11-S.13. For example, we can ask whether preparing a co-culture initially rich in cells of type z, as in (e), would generate a rapid changeover toward a population rich in cells of type x (horizontal gray arrow). We can investigate whether the population composition at (f) is a steady-state. We can also investigate whether preparing a co-culture with initial population composition at (g) would recapitulate the counterclockwise loop already traced out by the dataset in (c) or, instead, produce qualitatively different dynamics like the oscillating gray path. Such an oscillation could suggest the presence of a memory effect for reasons similar to those discussed in subsection 3.3 in the main manuscript.
Because our examples in Figure 3 from the main manuscript and Figure 5 happen to involve velocity fields and trajectories that resemble counterclockwise loops, it is helpful to inspect the examples in Figure 6. These velocity fields demonstrate that equations S.11-S.13 can produce a variety of features including curves and lines that look nothing like closed cycles. Additionally, the following section shows that EGT-based models and equations can accommodate the dynamics of biological populations in which cell-cell interactions are not an important effect. Thus, even when a set of replicator equations is 5/10 validated, careful scrutiny is necessary to determine whether the data could be adequately described using equations based on a simpler model of interaction-independent growth.

Parameter values
Replication rate coefficients substituted into equations 3.1 and 3.2 from the main manuscript and equations S.11-S.13 from the supplemental materials are organized below according to the figures they were used to generate. In this section, we also argue that careful scrutiny is needed to look for situations in which simple models of interaction-independent growth would suffice to describe an experimental system.  The parameter values in equations S.17 and S.18 happen to correspond to a prisoner's dilemma and a rough approximation of rock-paper-scissors, respectively. The names of these games are not important for the purposes of the tutorial.

Two-population examples
a differential equation describing population x expanding with growth coefficient -0.2/day (negative "expansion" actually corresponds to population decay), independent of population composition, and a differential equation describing population y expanding with growth coefficient 0.2/day, also independent of population composition.
This example shows that it is sometimes possible for EGT-based equations to fit population dynamics even when cell-cell interactions are not an important effect in an experimental system. In these cases, we cannot rely on a failed validation to alert us to the possibility of adequately interpreting data using simple models in which the cell-cell interactions so crucial to EGT are omitted. Figure 6(a) and (b) are similar, indicating that biological systems in which cell-cell interactions can be neglected can generate dynamics that qualitatively resemble the dynamics from some systems that do involve cell-cell interactions. A cursory inspection of a velocity field might not immediately reveal when an interactionindependent model of growth would be adequate. Taken together, this example shows that the dynamics of populations growing without cell-cell interactions (1) can pass the validation method described in subsection 3.2 and (2) can be described by phase portraits that look like phase portraits from systems in 7/10 which cell-cell interactions do contribute. Since validation and visual inspection of phase portraits might not alert us to the opportunity of using simple models in which cell-cell interactions are neglected, careful scrutiny of parameter values is required.