Cocktail, a Computer Program for Modelling Bacteriophage Infection Kinetics

Abstract

Cocktail is an easy-to-use computer program for mathematical modelling of bacteriophage (phage) infection kinetics in a chemostat. The infection of bacteria by phages results in complicated dynamic processes as both have the ability to multiply and change during the course of an infection. There is a need for a simple way to visualise these processes, not least due to the increased interest in phage therapy. Cocktail is completely self-contained and runs on a Windows 64-bit operating system. By changing the publicly available source code, the program can be developed in the directions that users see fit. Cocktail’s models consist of coupled differential equations that describe the infection of a bacterium in a vessel by one or two (interfering) phages. In the models, the bacterial population can be controlled by sixteen parameters, for example, through different growth rates, phage resistance, metabolically inactive cells or biofilm formation. The phages can be controlled by eight parameters each, such as different adsorption rates or latency periods. As the models in Cocktail describe the infection kinetics of phages in vitro, the program is primarily intended to generate hypotheses, but the results can however be indicative in the application of phage therapy.

1. Introduction

The ability of bacteriophages to kill bacteria has attracted increased interest in recent times. This is probably partly due to reports from clinics where they have been used for decades as bactericidal agents in the treatment of various infections, so-called phage therapy, but more recently also to a large number of studies of isolated phages (short for bacteriophages) which in laboratory experiments have been shown to be able to effectively kill bacteria in vitro, including antibiotic-resistant strains [1,2]. The increasing interest has also led to experimental treatments of severe infections, being carried out at other clinics than where they have traditionally been used for a long time [3,4]. The outcomes of clinical treatments are nevertheless difficult to predict, depending on many factors. Both bacteria and phages show an enormous variation between different clones and can change during the course of treatment, which together with largely unknown pharmacokinetics and pharmacodynamics results in treatments not being comparable and the effect of a certain treatment not being predictable. Although the variation between clinical experiments is very large, the combined results of them suggest that there is reason to continue studying how phage treatments against bacterial infections can be optimised.

A large number of mathematical models have been developed to study the dynamics between host organisms and their parasites, and bacteria–phages are no exception. Depending on the question, models use everything from evolutionary game theory (EGT), for example, on the evolution of the life cycle of phages, or networks of bacteria’s genetic changes after a bacteriophage infection (flux-balance analysis, FBA), but the most common bacteria–bacteriophage models are reaction kinetic models [5] (and references therein). These models essentially consist of coupled differential equations that describe the changes in the bacterial and bacteriophage population sizes over time in a vessel given a set of parameters. In their simplest form, the models describe the titre of susceptible, infected and resistant bacteria (classic SIR models) as well as the titre of the infecting bacteriophage. As mentioned above, bacteria and phages are self-replicating entities that can both change genetically, in the case of bacteria also phenotypically, through altered transcription [6,7]. Although simple models can provide valuable information, many models are more elaborate (for a review, see [8]). Several authors have contributed to the development of models that describe the dynamics of systems with, for example, more than one type of virulent bacteriophage, the formation of protective biofilms by bacteria or degradation of phages [9,10,11,12,13]. If models are also to be able to describe the dynamics in vitro, for example, during phage therapy of a human bacterial infection, it is also necessary to take into account the additional complexities that arise in interaction with human tissues and cells, for example, synergies between the immune system and phages, in the killing of bacteria or neutralisation of phages and phagocytosis by macrophages [14,15,16,17]. However, many of these in vitro interactions are still largely unknown in detail [18].

The purpose of the Cocktail program is to model the infection dynamics of one, or a combination of two, phage(s) infecting a bacterial species under varying relevant parameter settings, and in the latter case, to some extent, the interference between two phages infecting at the same time. The aim is to supply an easier way to carry out modelling in phage infection biology, e.g., for a better understanding of the complex dynamics during phage therapy. Although the program is easy to use, the most important parameters are included. The mathematical models are based on the basic models described, for example, by Levin et al. [10], Gill [13], Lenski [19], or Levin and Bull [20] and some of the additions discussed by Abedon [11]. It is thus possible in the program to also study the effect of metabolically inactive cells or biofilm formation as well as the decay of bacteria and phages, analogous to the action of an immune response (Figure 1).

As with all models, there has to be a balance between reality and generality. Hence, it is important to stress that some parameters, e.g., temperature, pH, release of nutrients from lysed bacteria or phages binding to cell debris, are not included in the program and others, e.g., modelling of the dynamics of an immune defence or cells in biofilm, are simplified. Therefore, the program should be seen as a tool for inducing hypotheses about the population dynamics of bacteria and phages during a phage infection and not exact predictions. This is especially important to stress regarding phage therapy experiments. The program does not consider the in vitro pharmacokinetics and pharmacodynamics of phages. Hypotheses will always need to be tested experimentally. Mathematical modelling of phage infection in a chemostat may however set boundaries to what can be expected while it can reflect the dynamics under ideal conditions (e.g., constant nutrient supply and agitation).

The following information describes the models and calculations in more detail. The basic condition is a bacterial population, growing in a vessel in a constant volume of nutrient, which can become infected by phages at varying titres and times. Although the volume is constant in such a chemostat, there could be an inflow and outflow of nutrients, and an outflow of bacteria and phages.

2. Materials and Methods

An overview of symbols used for all bacteria, phages and parameters can be found in

Table 1. Symbols and parameters.

		Start Values
Symbol	Description	Default	Allowed Range	Unit
Bacteria
S	Uninfected, susceptible bacteria	1 × 105	10−1 × 1012	CFU/mL
IA	Bacteria infected by phage A	-	-	|
IB	Bacteria infected by phage B	-	-	|
IAB	Bacteria infected by phages A and B	-	-	|
RA	Bacteria resistant to phage A	1 × 10−7	0–1 × 10−2	|
RB	Bacteria resistant to phage B	1 × 10−7	0–1 × 10−2	|
RAB	Bacteria resistant to phages A and B	1 × 10−14	0–1 × 10−6	|
RAIB	Bacteria resistant to A infected with B	-	-	|
RBIA	Bacteria resistant to B infected with A	-	-	|
Sr	Susceptible bacteria in a refuge	0	-	|
RrA	Bacteria resistant to A in a refuge	-	-	|
RrB	Bacteria resistant to B in a refuge	-	-	|
RrAB	Bacteria resistant to AB in a refuge	-	-	CFU/mL
Parameters
ψ	Growth rate of S	0.7	0–1.5	/h
K	Monod constant	5.0	0.01–100	µg/mL *
ε	Resource for division of one bacterium	2 × 10−6	1 × 10−8–1 × 10−4	µg/cell *
γ	Bacterial decay rate	0	0–1	/h
µA	Mutation rate for resistance against A	1 × 10−7	0–1 × 10−2	/cell div.
µB	Mutation rate for resistance against B	1 × 10−7	0–1 × 10−2	/cell div.
${\psi }_{R_A}$	Growth rate of RA	0.7	0–1.5	/h
${\psi }_{R_B}$	Growth rate of RB	0.7	0–1.5	/h
${\psi }_{R_{AB}}$	Growth rate of RAB	0.7	0–1.5	/h
σ	Rate of bacteria into refuge	0	0–0.01	/min
ρ	Rate of bacteria out from refuge	0	0–0.01	/min
C0	Available resources from start	100	0–1000	µg/mL *
C	Resources flowing in from a reservoir	100	0–1000	µg/mL *
ω	Flow rate	0.2	0–100	/h
Phages
Parameters
A	Titre of phage A	1 × 108	0–1 × 1013	PFU/mL
B	Titre of phage B	1 × 108	0–1 × 1013	PFU/mL
δA	Adsorption rate of A	1 × 10−10	1 × 10−14–1 × 10−7	mL/min
δB	Adsorption rate of B	1 × 10−10	1 × 10−14–1 × 10−7	mL/min
lA	Latent period of A	30	1–60	min
lB	Latent period of B	20	1–60	min
βA	Burst size of A	100	0–1000	PFU/cell
βB	Burst size of B	100	0–1000	PFU/cell
φA	Decay rate of phage A	0	0–1	/h
φB	Decay rate of phage B	0	0–1	/h

* The symbol for the micro prefix, “µ”, is denoted by “u” in the program user interface.

. In the program, values can in general be given with three significant digits. If the input should be an integer, it can be given either as that or in scientific notation, e.g., 1,000,000 or 1.0 × 10⁶. Real numbers should be given either in decimal or scientific format with a point as the decimal separator. If the wrong format is used, values are auto corrected in most cases. Hovering with the mouse pointer above a box displays a hint on which values can be given. The tab key can be used for jumping between boxes and it is also possible to use the up and down arrows in some boxes to increase or decrease values in fixed steps. Please note that some parameter values should be given per hour and others per minute, as shown by the default values (Table 1). Additionally, note that many parameters are represented with their average values despite in reality being distributed in time or size, e.g., the latent period will vary from cell to cell and not all cells will generate the same burst size (produce exactly the same number of phages).

2.1. Bacteria

For most bacteria, the rate of growth mainly depends on the concentration of nutrients (physical factors, e.g., temperature and the presence of various gases are of course also important). In a closed system, while nutrients are consumed by the bacteria, their concentration decreases and growth slows. Seen over time, in such cases the bacterial growth becomes a logistic function of the concentration of nutrients. One of the most widely used relationships between growth rate and nutrient concentration was formulated by Monod [21]; the growth of bacteria depends on three parameters, µ = µ_max × s/(K_s + s) where the growth rate, µ, is a function of the maximum growth rate, µ_max, a constant, K, and the concentration of nutrients, s. This is also the basic growth function in the Cocktail program where most parameter symbols and default parameter settings are taken from Lenski [19] (Table 1). In the program, the max growth rate per hour is denoted as ψ, a real number between 0 and 1.5 and the Monod half-saturation constant, K in µg/mL, is a real number between 0.01 and 100. The reason for allowing this large span is the observed variation among strains and experimental conditions when assessing the half-saturation constant, but the values are usually between 0.1 and 10 µg/mL [22]. The concentration of nutrients varies in a chemostat while the nutrients are consumed by bacteria and where there is an inflow and outflow. It is denoted C in the program and the concentration of nutrients over time equation is described in more detail below. While mutated bacteria may suffer from reduced fitness, the growth rate of bacteria that have become resistant to either one or both phages can be entered separately. Bacteria can become resistant to phage infection by mutation but only at cell division when new cells are produced. The mutation rates should also be written as decimal numbers or in scientific notation. Note that bacteria becoming resistant to both phages, A and B, occurs at a rate being the product of the rate of mutation to resistance against A and B, respectively, and does not need to be specified. It is also possible that the starting population of bacteria may contain resistant bacteria to either one or both phages. These frequencies are also entered either as decimal numbers or in scientific notation as above.

Bacteria can also decay from natural causes, and not just die from a phage infection. In an in vitro situation, this would be, e.g., from neutralisation by antibodies or phagocytosis. This results in an exponential decay N_t = N₀ × e^−γt where N₀ is the number of bacteria at time 0 and N_t at time t and γ is the decay rate constant. The decay is calculated as N_t+1 = N_t − (N_t × γ) in the equations below while this equals N₀ × e^−γt. γ is quite small for bacteria growing in a chemostat, half of the population has decayed at t = ln(2)/γ, and values of γ is generally around 0.02 per hour in nature [23,24]. The decay rate per hour, γ, should be given as a real number between 0 and 1. Note that new bacteria resulting from cell division each generation is excluded from decay.

2.2. Resources

The addition of nutrition is necessary for the growth of bacteria. The concentration of nutrients is regulated by the flow into and out of the system at a rate of ω turnovers/hour. The starting concentration, C₀ in µg/mL, can be set independently from the concentration of nutrients continuously flowing into the system, C in µg/mL. The flow in and out of the system causes dynamic changes in the system while it is not just nutrients flowing out but also bacteria and phages. The conversion efficiency, the resources used by one dividing cell, is denoted as ε and given in µg/cell as a decimal number or in scientific notation, e.g., 1.0 × 10⁻⁶. The equation for dC/dt is shown below in the following text.

2.3. Phages

Two virulent phages, A and B, can infect the bacterial population in different titres, adsorption rates, and at different times while it is possible to let the bacterial population grow for some time and add phages at three different time points. As with the bacteria, phages can decay and be washed out of the system. Phages inactivated by binding to lysed cell debris cannot be entered separately but must be included in the decay. Phage populations grow by two additional parameters. When the latent period comes to an end, infected bacteria burst open producing the burst size number of phages. Note that a fast phage of one type, e.g., A, infecting a bacterium already infected with a slower phage, e.g., B, will produce only type A phages if its latent period is shorter than the remaining time of phage B’s latent period. If only one phage is chosen to infect, the added titre of the other phage should be set to 0. Additionally, if one wants to study a nonproductive infection, the burst size should be set to zero.

2.4. Model Settings

A phage infection starts with the adsorption of the phage to receptors on the bacterial surface. The adsorption can be set to happen by different mathematical models described in the following text. In short, phages can adsorb one by one per unit time as described in the primary adsorption “Standard” model below or several phages at once per unit time according to the Poisson probability of infection. The “Poisson” setting hence allows for multiple adsorption proportional to the multiplicity of infection. It is also possible to adjust both models to allow adsorption to uninfected and non-resistant cells only or to all susceptible and non-resistant cells irrespective if they have already been infected. These secondary adsorption alternatives are chosen with the settings “Uninfected” or “Susceptible”. An overview over which types of cells that get adsorbed and infected by which phages is given in

Table 2. Adsorbed A and B phages per unit time under the primary and secondary adsorption settings.

Primary Adsorption Setting	Standard		Poisson
Secondary Adsorption Setting	Uninfected	Susceptible	Uninfected	Susceptible
	Phages adsorb one at a time to uninfected non-resistant cells	Phages adsorb one at a time to non-resistant cells	A number of phages adsorb according to a Poisson probability with lambda = MOI to uninfected non-resistant cells	A number of phages adsorb according to a Poisson probability with lambda = MOI to non-resistant cells
Bacteria	Conceivably adsorbing phages
S = Susceptible	A or B	A or B	A and B	A and B
IA = Infected by A	B	A or B	B	A and B
IB = Infected by B	A	A or B	A	A and B
IAB = Infected by A and B	-	A or B	-	A and B
RA = Resistant to infections by A	B	B	B	B
RB = Resistant to infections by B	A	A	A	A
RAB = Resistant to infections by A and B	-	-	-	-
RAIB = Resistant to infections by A infected with B	-	B	-	B
RBIA = Resistant to infections by B infected with A	-	A	-	A
Sr = Susceptible planktonic bacteria in a refuge	-	A or B No infection	A and B No infection	A and B No infection
RrA = Planktonic bacteria resistant to A in a refuge	-	B No infection	B No infection	B No infection
RrB = Planktonic bacteria resistant to B in a refuge	-	A No infection	A No infection	A No infection
RrAB = Planktonic bacteria resistant to AB in a refuge	-	-	-	-
Sr = Susceptible bacteria in a LIFO refuge	-	-	-	-
RrA = Bacteria resistant to A in a LIFO refuge	-	-	-	-
RrB = Bacteria resistant to B in a LIFO refuge	-	-	-	-
RrAB = Bacteria resistant to AB in a LIFO refuge	-	-	-	-

. Details and the mathematical background to all models is given in the text that follows.

2.4.1. Primary Adsorption: Standard Model

Bacteria can divide and the population grow, but bacteria can also decay, be completely resistant or mutate to resistance against phage infection (Section 2.4.3), hide in a refuge (Section 2.4.4) as well as become infected by phages and lyse from the infection. In the Cocktail program, this results in the possibility of thirteen types of bacteria being present in the system at the same time (Table 1). As mentioned, phages A and B can infect at different times and by different adsorption rates, and after varying latent periods lyse infected bacteria resulting in a burst of phages of different size. Bacteria and phages can both flow out of the system. Starting with the concentration of nutrients, this results in the following basic equations:

$$\frac{dC}{dt}=\left(C_0-C\right)\omega -\frac{\epsilon N\psi C}{K+C}$$

(1)

The first Equation (1) describes the nutrient content in the system over time where C is the concentration of nutrient in µg/mL (C₀ is the start concentration), ω is the flow rate in turnovers/h, ε is the resource consumption by one bacterium in µg/cell and ψ is the specific growth rate of bacteria per hour. N stands for the S and R types of bacteria (S, R_A, R_B, R_AB) as the infected bacteria are supposed to cease both to consume nutrients and to divide. K is the Monod half-saturation constant in µg/mL (the concentration of nutrients resulting in half the maximum growth of the bacteria). For simplicity, ε and K are the same for all dividing bacteria, but ψ can vary and be different for phage-resistant bacteria.

$$\frac{dS}{dt}=\frac{S{\psi }_{S}C}{K+C}-\frac{({\mu }_A+{\mu }_B+{\mu }_{AB)}S{\psi }_{S}C}{K+C}+\rho S_r-\sigma S-\mathsf{\gamma }S-{\delta }_{A}SA-{\delta }_{B}SB-\omega S$$

(2)

This second Equation (2) describes the growth and losses of susceptible bacteria. The losses over time are due to bacteria mutating to become resistant against phage A or B or both at a specific rate of μ per cell division (µ_AB is calculated as µ_A × µ_B). This means that it is only divided bacteria that mutate. Bacteria can also move into a refuge population, S_r, at a rate of σ where they may or may not be adsorbed by phages, and move out of the refuge at a rate of ρ. Bacteria can also decay or be neutralised at a rate γ, becoming infected by phages at the adsorption rate δ, and finally be washed out of the system at a rate of ω. The equations for phage-resistant bacteria become:

$$\frac{d{R}_A}{dt}=\frac{R_A{\psi }_{R_A}C}{K+C}+{\mu }_{A}S+\rho R_{rA}-\sigma R_A-\mathsf{\gamma }{R}_A-{\mu }_B{R}_A-{\delta }_B{R}_{A}B-\omega R_A$$

(3)

$$\frac{d{R}_B}{dt}=\frac{R_B{\psi }_{R_B}C}{K+C}+{\mu }_{B}S+\rho R_{rB}-\sigma R_B-\mathsf{\gamma }{R}_B-{\mu }_A{R}_B-{\delta }_A{R}_{B}A-\omega R_B$$

(4)

$$\frac{d{R}_{AB}}{dt}=\frac{R_{AB}{\psi }_{R_{AB}}C}{K+C}+{\mu }_{AB}S+\rho R_{rAB}-\sigma R_{AB}-\mathsf{\gamma }{R}_{AB}+{\mu }_A{R}_B+{\mu }_B{R}_A-\omega R_{AB}$$

(5)

Phage resistance is intended to result in complete blocking of adsorption by the phage. Hence, a bacterium resistant to phage A can mutate to become resistant to phage B as well but also be infected by phage B, and vice versa. Needless to say, a bacterium resistant to both phages (Equation (5)) cannot become infected at all.

Infected bacteria (Equations (6)–(8)) are thought not to consume resources or divide and not to be part of the refuge population of cells. Cells in the refuge (see Equations (26) and (27) below) is simulating the presence of either metabolically inactive cells or cells in biofilm and inhibited phage propagation. Phages can however adsorb to different classes of bacteria depending on the secondary adsorption mode setting. In the “Uninfected” setting, phages are adsorbing to uninfected bacteria only, which is common in basic mathematical models, where phage A adsorbs only to S, I_B, and R_B. In the “Susceptible” setting, on the other hand, phages are allowed to adsorb to already-infected bacteria as well, i.e., A adsorbs to S, I_A, I_B, I_AB, R_B and R_BIA (referred to as secondary adsorption [11]). In addition to this, A can also adsorb to the S_r and R_rB cells in the refuge if the “Standard” primary adsorption model is applied, the “Planktonic” mode is set and the rate of cells in and out of the refuge are given values. Phage B adsorbs to S_r and R_rA as well with this setting. Cells in the refuge will however not be infected, only adsorbed by phages.

$$\frac{d{I}_A}{dt}={\delta }_{A}SA-\mathsf{\gamma }{I}_A-{\delta }_B{I}_{A}B-I_{A\left(t-l_A\right)}-\omega I_A$$

(6)

$$\frac{d{I}_B}{dt}={\delta }_{B}SB-\mathsf{\gamma }{I}_B-{\delta }_A{I}_{B}A-I_{B\left(t-l_B\right)}-\omega I_B$$

(7)

$$\frac{d{I}_{AB}}{dt}={\delta }_A{I}_{B}A+{\delta }_B{I}_{A}B-\mathsf{\gamma }{I}_{AB}-I_{AB\left(t-l_B\right)}-I_{AB\left(t-l_A\right)}-\omega I_{AB}$$

(8)

The expressions $I_{A\left(t-l_A\right)}$, $I_{B\left(t-l_B\right)}$, $I_{AB\left(t-l_B\right)}$ and $I_{AB\left(t-l_A\right)}$, part of the delayed differential Equations (6)–(8), are all describing the loss of infected bacteria due to lysis, after the latency time l_A or l_B. A bacterium simultaneously infected by two phages will lyse at time l_A if the latency of phage A is shorter than the latency time for phage B, l_B. This results in that $I_{AB\left(t-l_B\right)}$ becomes zero for a bacterium when $I_{AB\left(t-l_A\right)}$ becomes positive. One of $I_{AB\left(t-l_A\right)}$ and $I_{AB\left(t-l_B\right)}$ subsequently always becomes zero. A bacterium infected with phage B at time t_B and superinfected with phage A at time t_A will lyse and produce phages of type A only if t_A + l_A is shorter than the remaining time to lysis caused by phage B, i.e., t_A + l_A < t_A + l_A − t_B + l_B = t_A < t_A − t_B + l_B. This means that in an infection with both phages, A and B will interfere with each other and not produce the number of phages expected from single and separate infections by A or B. Other interference between phages, e.g., by actively impeding the other phage’s transcription or replication, is not considered.

Finally, bacteria that are resistant to infections by one phage can become infected with the other phage (Equations (9) and (10)).

$$\frac{d{R}_{BIA}}{dt}={\delta }_A{R}_{B}A-\mathsf{\gamma }{R}_{BIA}-{\delta }_A{R}_{BIA\left(t-l_A\right)}-\omega R_{BIA}$$

(9)

$$\frac{d{R}_{AIB}}{dt}={\delta }_B{R}_{A}B-\mathsf{\gamma }{R}_{AIB}-{\delta }_B{R}_{AIB\left(t-l_B\right)}-\omega R_{AIB}$$

(10)

Titres of phages A and B can grow through the release of phages from all types of infected bacteria (13 and 14). After the phage specific latency periods mentioned above, each bacterial cell gives rise to the number of phages equal to the phage’s burst size, β. Phages will be lost by adsorption to the bacteria mentioned above, and adsorbed phages representing phages bound to bacteria, is denoted P_A and P_B, respectively (11 and 12). Phages can also decompose at a rate of φ, and be washed out at a rate of ω. With the no secondary adsorption “Uninfected” setting:

$$P_A={\delta }_{A}A\left(S+I_B+R_B\right) \mathrm{and} P_B={\delta }_{B}B\left(S+I_A+R_A\right)$$

(11)

Additionally, with the “Susceptible” setting:

$$\begin{gathered} P_A={\delta }_{A}A\left(S+I_B+R_B+I_A+I_{AB}+R_{BIA}\right) \\ \mathrm{and} \\ {P}_B={\delta }_{B}B\left(S+I_A+R_A+I_B+I_{AB}+R_{AIB}\right) \end{gathered}$$

(12)

The change in phage titres will hence be:

$$\frac{dA}{dt}={\beta }_A\left(I_{A\left(t-l_A\right)}+I_{AB\left(t-l_A\right)}+R_{BIA\left(t-l_A\right)}\right)-P_A-{\phi }_{A}A-\omega A$$

(13)

$$\frac{dB}{dt}={\beta }_B\left(I_{B\left(t-l_B\right)}+I_{AB\left(t-l_B\right)}+R_{AIB\left(t-l_B\right)}\right)-P_B-{\phi }_{B}B-\omega B$$

(14)

2.4.2. Secondary Adsorption: Poisson

Only uninfected bacteria can be infected by phages in the “Standard” model with the “Secondary adsorption Uninfected” setting (${\delta }_{A}SA$ in the equations). This is in many cases a good approximation but neglects that several phages may adsorb to a single bacterial cell, assuming that the number of cell receptors is not limited, i.e., multiple adsorption [11]. Phages adsorbed to cells will follow Poisson probabilities. While more phages per bacterium can infect, this is referred to as MOI_actual in contrast to MOI_input [11,25]. Bound phages will hence be:

$$P_{t+1}=\left(1-e^{-\delta M_t}\right)P_t$$

(15)

Here, P is the titre of phages and M is the sum of all bacteria that can be adsorbed by a particular phage. These are denoted M_A for phage A and are bacteria S, I_B, and R_B, in the “Uninfected” adsorption mode and S, I_A, I_B, I_AB, R_B and R_BIA with the “Susceptible” setting, i.e., the same sets of bacteria as in the “Standard” primary adsorption model. M_B is accordingly bacteria S, I_A, and R_A, and S, I_B, I_A, I_AB, R_A and R_AIB, respectively. While more phages than one may adsorb to a cell, infected bacteria will equal:

$$I_{t+1}=(1-e^{-\frac{P_t}{M_t}})M_t$$

(16)

Taking Equation (16) into account, the equations describing the change in uninfected bacterial titres (17–19) will have to change to:

$$\frac{dS}{dt}=\frac{S{\psi }_{S}C}{K+C}-\frac{({\mu }_A+{\mu }_B+{\mu }_{AB)}S{\psi }_{S}C}{K+C}+\rho S_r-\sigma S-\mathsf{\gamma }S-((1-e^{-\frac{P_A}{M_A}})M_A)-((1-e^{-\frac{P_B}{M_B}})M_B)-\omega S$$

(17)

The equations for resistant bacteria will change accordingly:

$$\frac{d{R}_A}{dt}=\frac{R_A{\psi }_{R_A}C}{K+C}+{\mu }_{A}S+\rho R_{rA}-\sigma R_A-\mathsf{\gamma }{R}_A-{\mu }_B{R}_A-((1-e^{-\frac{P_B}{R_A}})R_A)-\omega R_A$$

(18)

$$\frac{d{R}_B}{dt}=\frac{R_B{\psi }_{R_B}C}{K+C}+{\mu }_{B}S+\rho R_{rB}-\sigma R_B-\mathsf{\gamma }{R}_B-{\mu }_A{R}_B-((1-e^{-\frac{P_A}{R_B}})R_B)-\omega R_B$$

(19)

$$\frac{d{R}_{AB}}{dt}=\frac{R_{AB}{\psi }_{R_{AB}}C}{K+C}+{\mu }_{AB}S+\rho R_{rAB}-\sigma R_{AB}-\mathsf{\gamma }{R}_{AB}+{\mu }_A{R}_B+{\mu }_B{R}_A-\omega R_{AB}$$

(20)

If the infection is by the primary adsorption option “Poisson” and bacteria are chosen to be in the planktonic refuge and entered at a certain rate (see Equation (26) below), phages are additionally also adsorbing to the bacteria in the refuge. Phage A will additionally adsorb to S_r and R_rB cells and phage B to S_r and R_rA cells. However, this is not the case if the last-in-first-out (“LIFO”) type of refuge cells is selected (27). Bound phages, A and B, are again denoted P_A and P_B, respectively. In the Poisson mode, this results in Equations (21)–(25) for infected bacteria:

$$\frac{d{I}_A}{dt}=((1-e^{-\frac{P_A}{M_A}})M_A)-\mathsf{\gamma }{I}_A-{\mu }_B{I}_A-((1-e^{-\frac{P_A}{I_B}})I_B)-I_{A\left(t-l_A\right)}-\omega I_A$$

(21)

$$\frac{d{I}_B}{dt}=((1-e^{-\frac{P_B}{M_B}})M_B)-\mathsf{\gamma }{I}_B-{\mu }_A{I}_B-((1-e^{-\frac{P_B}{I_A}})I_A)-I_{B\left(t-l_B\right)}-\omega I_B$$

(22)

$$\frac{d{I}_{AB}}{dt}=((1-e^{-\frac{P_A}{I_B}})I_B)+((1-e^{-\frac{P_B}{I_A}})I_A)-\mathsf{\gamma }{I}_{AB}-I_{AB\left(t-l_B\right)}-I_{AB\left(t-l_A\right)}-\omega I_{AB}$$

(23)

$$\frac{d{R}_{BIA}}{dt}=((1-e^{-\frac{P_A}{R_B}})R_B)-\mathsf{\gamma }{R}_{BIA}-R_{BIA\left(t-l_A\right)}-\omega R_{BIA}$$

(24)

$$\frac{d{R}_{AIB}}{dt}=((1-e^{-\frac{P_B}{R_A}})R_A)-\mathsf{\gamma }{R}_{AIB}-R_{AIB\left(t-l_B\right)}-\omega R_{AIB}$$

(25)

The expression of phages A and B lost by adsorption to bacteria is the same as in the “Standard” model, described in Equations (13) and (14), but bound phages P_A and P_B is calculated differently as shown in Equation (15) above.

It should be pointed out that the difference between the “Standard” mode of infection and the “Poisson” mode is obviously small at a MOI around 1. It is only when there are phages in excess, a probability of more than one phage infecting a bacterium, that a difference may be observed as a more rapid loss of adsorbed phages.

2.4.3. Resistance Mutation

Dividing bacteria can mutate at a rate of µ, and the mutant frequencies in the population can be calculated in two different ways. In the “Deterministic” mode, each class of newly divided bacteria contains N × µ mutants, where N is the number of newly divided bacteria. With the “Stochastic” alternative, mutations are introduced by random sampling from a Poisson distribution having a mean of N × µ = λ if λ ≤ 10 or from sampling a normal distribution if λ > 10. The normal distribution, (λ;√ λ), is generated by the Box-Müller algorithm and all random numbers are generated by a Mersenne Twister algorithm. This results in good, but somewhat slow, generation of pseudo random numbers but this does not have a great impact on the overall program performance. When stochastic mutation is active, results will of course vary from run to run. With small numbers of divided bacteria per millilitre, for example, 10⁵ bacteria, and a mutation rate of 10⁻⁷, the resulting frequency of mutants is bound to be very low. There would be only 0.01 mutant bacteria in the population and these would be eliminated if the option of rounding off values below one is activated. In the program, resistance is modelled as affecting the adsorption and regarded as complete. Therefore, resistant bacteria do not adsorb any phages.

2.4.4. Refuge Cells

The refuge population is simulating the existence of metabolically inactive cells, with the “Planktonic” setting, or cells forming biofilm with the last-in-first-out (“LIFO”) setting. The “LIFO” setting means that the last cells that entered the refuge are reintroduced to the normal pool of cells followed by the next to last cells and so forth. The rate of cells moving into the refuge can be set to at most 0.01, which means that 1% of the current population will enter the refuge per minute. Another limitation is that there has to be more than 10 cells outside of the refuge. In such a case, only 0.1 cells enter the refuge. If only whole cells should be allowed to enter, the round off <1 option should be chosen. Both boxes, the rate in/min and rate out/min must be given a value in order to activate the refuge cells models. While in the refuge, these cells are not dividing and cannot produce phages or mutate and become resistant which they can become upon reintroduction to the normal pool of cells. Infected bacteria (I_A, I_B, I_AB, R_AIB, R_BIA) are not part of either refuge population, as these cells will lyse in any event. In the “Planktonic” mode, cells can be flushed out of the system or decay whereas in “LIFO” mode cells are thought to be metabolically inactive and sessile until reintroduced. Planktonic refuge cells:

$$\begin{gathered} \frac{d{S}_r}{dt}=\sigma S-\rho S_r-\mathsf{\gamma }{S}_r-\omega S_r ; \frac{d{R}_{rA}}{dt}=\sigma R_A-\rho R_{rA}-\mathsf{\gamma }{R}_{rA}-\omega R_{rA} \\ \frac{d{R}_{rB}}{dt}=\sigma R_B-\rho R_{rB}-\mathsf{\gamma }{R}_{rB}-\omega R_{rB} ; \frac{d{R}_{rAB}}{dt}=\sigma R_{AB}-\rho R_{rAB}-\mathsf{\gamma }{R}_{rAB}-\omega R_{rAB} \end{gathered}$$

(26)

Cells in the last-in-first-out (LIFO) mode:

$$\begin{gathered} \frac{d{S}_r}{dt}=\sigma S-\rho S_r ; \frac{d{R}_{rA}}{dt}=\sigma R_A-\rho R_{rA} \\ \frac{d{R}_{rB}}{dt}=\sigma R_B-\rho R_{rB} ; \frac{d{R}_{rAB}}{dt}=\sigma R_{AB}-\rho R_{rAB} \end{gathered}$$

(27)

2.4.5. Time Step Size

Calculations of differential equations in Cocktail are done using Euler’s method, taking the input values as the initial values for calculating the values numerically after the chosen time interval, set either as one minute or as a 30-, 15- or 5-second time step size. At large time intervals, other methods for solving differential equations result in smaller errors, but the differences between methods become smaller and smaller as the step size (time interval) decreases. Hence, running the program with a step size of one minute results in a larger discretisation error than with a five-second step size, but is considerably faster. On the other hand, the program speed depends mainly on the length of the phages’ latent periods and the rates of bacteria into and out of refuge populations. These are stored on arrays in dynamic memory that need to be recalculated in each step, which will slow down the program performance if a long running time is set. However, there is virtually no time difference between short and long time step sizes if a short running time and no refuge cells are chosen. If accuracy and a small discretisation error is preferred, step size should be set to five-second intervals.

3. Results

Examples of Cocktail outputs can be studied by running the example data files provided as Supplementary Materials. The file Bohannan_Lenski_1997_Fig 3B.ctl contains input parameter values and a comparison to a chemostat experiment where the authors found that the titre of Escherichia coli bacteria and an infecting T4 phage can oscillate over time [26]. The phage titre decreases over time when there are very few bacteria to infect which in turn results in a higher titre of bacteria and so forth (Figure 2A). Coexistence can theoretically be shown to occur at higher bacterial titres as well. Running the parameter settings in the file Lenski_1988_Fig_2a.ctl from Lenski [19] results in oscillations leading to stable coexistence of bacteria in a titre of about 10⁷ and phage titres being around 10⁹ (Figure 2B). However, a fourfold increase in the concentration of nutrients, from 25 to 100 µg/mL results in increasingly large oscillations, but as in the first case, bacteria never become extinct (Figure 2C).

It is also possible to analyse more complex problems and formulate hypotheses that later can be evaluated experimentally. The last example describes two phages in a cocktail with different latent periods, added in the same titres, and at the same time to the bacteria. Most parameters were entered with their default values (Table 1). Parameters set to different values were the bacteria’s starting titre, 1 × 10⁸, and the program was executed with the standard model where phages only adsorb to uninfected cells, mutations are set to be deterministic and no cells were allowed to be resistant to either phage A or B or both, metabolically inactive or form biofilm. The log₁₀ option was selected for the output. While log₁₀(0) = −∞, a titre of 0 is represented as −16 (log₁₀ of 10⁻¹⁶). The data can be retrieved by running the file Fig_2D.ctl. The result showed that bacteria resistant to phage A will disappear from the system within an hour (Figure 2D). Phage B has a shorter latent period, and has outcompeted phage A, by infecting most of the susceptible bacteria. The mutation rate of bacteria becoming resistant to both phages was set to 10⁻⁷ × 10⁻⁷ (the product of mutation rate for resistance to A and B, respectively) which resulted in a low titre of double resistant cells, as the bacteria had a good supply of nutrients and were able to divide, but the titre of such cells will grow to about a thousand cells/mL in 48 h. Allowing resistant cells from the beginning, in the start population, results in higher titres of such cells.

4. Technical Information

Cocktail runs on Windows 64-bit systems. The program interface (Figure 1) is in English, but some instructions may turn up in the language set on your computer when Windows DLLs are called (e.g., file dialogs). There is no support for other languages in Cocktail. The program is developed in Object Pascal from the Free Pascal Team (Free Pascal: A 32-, 64- and 16-bit professional Pascal compiler. Version 3.2.0. URL https://www.freepascal.org. RRID:SCR_014360), accessed on 28 September 2022, using the Lazarus IDE and libraries developed by the Lazarus Team (Lazarus: The professional Free Pascal RAD IDE. Version 2.0.10. URL http://www.lazarus-ide.org. RRID:SCR_014362, accessed on 28 September 2022). The IDE, compiler and program libraries can be downloaded from: https://www.lazarus-ide.org/index.php?page=downloads, accessed on 28 September 2022. Source code files (cocktailunit1.pas, cocktailunit2.pas, cocktailunit1.lfm, cocktailunit2.lfm, Cocktail.lps, Cocktail1.lpr, Cocktail1.res, Cocktail.dbg, globalvariables.pas, Cocktail.ico and the Lazarus project information file Cocktail.lpi) as well as updates will be available from GitHub at https://github.com/ASNilsson/Cocktail-phage-infection-kinetics, accessed on 28 September 2022.

The Cocktail program and source code files are distributed under the license Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International License. In short, this means that it is free for everyone to use, to modify the source code, build upon the program or code, and free to distribute in any medium. However, you must give appropriate credit and a link to the license. If changes were made to the program or code, these must be specified, and distribution of modifications must be under the same license. It is not allowed for anyone to use any part of the program or code for commercial purposes. A short description of the license can be found at: https://creativecommons.org/licenses/by-nc-sa/4.0/, accessed on 28 September 2022. The license and program version number can be found by double clicking anywhere in the Cocktail parameter settings window.

The results of the program can be saved as charts, in PNG or SVG graphics file formats, and/or as a .ctl data file. The items in the .ctl file constitute the complete settings for running the program. The file format is a plain .txt file, but note that the format is fixed as in the example .ctl file. Moving items to another position in the file will inevitably result in a file error. A comma (,) is often used as a delimiter when more than one item is to be found on a line. The advantage of this is that the file can contain a short label of each of the items which makes reading and editing a file much easier. The disadvantage is that omitting a comma, or using another delimiter, will result in a file error. Selected output parameters should however be surrounded by at least one blank in the list following the label “Output parameters:”, e.g., 1 11 14 (note the blank after the last number). A .ctl file can easily be created by running the program and saving the result by clicking on the “Save” button at the bottom of the result window. Editing such a .ctl file that has been shown to work as a template, and saving it under a new name, is a good idea. Double clicking on a .ctl file will open a new instance of the program provided that a link to the program has been established in the Windows “How do you want to open this file?” dialog by marking the Cocktail program and checking the “Always use this app to open .ctl files” box.