Model: Methodology

Overall approach

In order to understand the interactions between phage and bacteria, and build a model that would incorporate both temperate and lytic phages, all models mentioned in the background had been built and simulated separately (Campbell, 1961; Levin, Stewart and Chao, 1977; Qiu, 2007). Although intermediate results are not presented here, one can find source files in our repository on GitHub. Simulations of the final models are presented in Results. All these models are deterministic, meaning that the outcome of the model does not change as long as the input values are kept the same. This does not impose issues from the mathematical point of view, however, physically it makes less sense as the values approach low numbers. Thus, these models fail to accurately represent the true experimental conditions for low values of XS, XI and P (Smith and Trevino, 2009).

The chemostat was used in all previously models for the continuous culture. This is due to the following reasons: (1) it is considered as the best apparatus that allows imitation of nature in laboratory for population studies (Hsu, Hubbell and Waltman, 1977), (2) it allows simple manipulations and control over experimental parameters and, (3) chemostat models usually show consistent results with experimental data (Hansen and Hubbell, 1980).

Simulations

Python programming language was used for scripting and simulations. In particular, an open-source SciPy (0.19.0) software was used for general numerical manipulation of values and matplotlib.pyplot for plotting in Python. Although the most widely used package for integration of a system of ODDs odeint is included into SciPy by default, its functionality was not enough in order to solve systems of DDEs. Therefore, a side package called PyDDE (0.2.2) was used to solve these systems of equations.

In contrast to ODEs, ODDs allow to break out from the continuous calculations and take into account the value of the system from the past. This is particularly relevant for modelling of phage infections due to the latency time. Therefore, the number of phages that are released at time t depends on the number of bacteria that got infected at t-T, where T is the time latency.

Source code

All code used for the simulations is open-source and stored in version control repository Github (https://github.com/csynbiosys/iGEM-Phage-Modelling/). For convenience, every plot generated in this work in the title contains a link to a hash (SHA) of a particular commit on Github. This allows every result presented here to be easily reproduced using Python - one of the easiest programming languages available.

Calculation of the r-value

The average concentration was calculated by by adding up all concentrations above 10-15 particles/ml and dividing this by the number of concentrations points taken. This was done for concentrations of lytically infected and susceptible bacteria. The value r was then calculated as a ratio of XI to XS.

Initial values

The majority of the initial values for the simulations were taken from Levin et al. (1977) and are listed in Table 1. The carrying capacity was calculated as for Escherichia coli in LB growth media by multiplying optical density at the carrying capacity (OD600 ≈ 7) by the number of bacteria per 1 OD600 (109 cell/mL), resulting in 7.0·109 cells. The maximum growth rate of E.coli in presence of 0.2% L-arabinose was 0.568 per hour (Aidelberg et al., 2014). The half-saturation constant (Km) for L-arabinose was taken as five times higher than that of glucose (Meyenburg, 1971) i.e. 20 μg/ml.

Table 1. Parameters used in models of populations of phage and bacteria.

Sensitivity analysis

We performed a local sensitivity analysis on all 20 parameters used in the model by calculating the elasticity, which shows the percentage change in the output variable divided by the percentage change in the parameter. Two main output variables were tracked: the extinction time of XS and r-value. As the main objective of the sensitivity analysis was to quantify the effect of each variable on the key outputs, we simply varied the parameters from the minimum value at 50% of the nominal value and for the maximum value at 150%. The standard formulae for the calculation of elasticity is as follows:

, where F – is the output of the objective function when the parameter is either nominal (n) or alternative (a) and similarly P – is the parameter. In our case the change percent change in the parameter is already known, therefore, the formulae above is simplified to:

, where F – is the output of the objective function when the parameter is either 150% or 50% of the nominal value.

The elasticity is a good and simple measure of sensitivity, however it gives reliable information only if the two variables of interest have a linear relationship. Therefore, we calculated Pearson Correlation Coefficient for every parameter of the model along the span between 50% and 150% of the nominal.