iGEMwiki toolssearchteamsdoucheno1

Editing Team:IISc-Bangalore/Model

Modelling

Overview

Bacteria–phage coevolution, the reciprocal evolution between bacterial hosts and the phages that infect them, is an important part of this year's project. There is growing evidence from both laboratory and natural populations that coevolution can maintain phenotypic and genetic diversity, increase the rate of bacterial and phage evolution and divergence, affect community structure, and shape the evolution of ecologically relevant bacterial traits.^[6] Also, phage dosage has been an elusive and non-trivial concept, proving to be a considerable challenge in the field of therapeutics. Hence we decided to have a go at it by employing delayed ODE-based methods. Additionally, using a stochastic approach, we modeled induced mutations in bacteriophage T4 to determine feasibility of obtaining mutants capable of binding to resistant bacteria.

Two models were devised in the course of the project.

Numerically analysing the phage-bacteria interaction using modified predator-prey equations.
Calculating the expected number of a particular genome sequence arising from the T4 phage mutations.

Coevolution

Due to the self-replicating nature of phages, some peculiarities arise in the use of phages as an antimicrobial drug. Two modes of therapy are possible with phages, active and passive. Active therapy takes place when the virus population increases along with the bacterial population, eventually overwhelming the bacterial population and drastically reducing bacterial concentration. Passive therapy is when the initial viral population is sufficiently large that the bacterial population is wiped out by primary infection alone, and displays kinetics somewhat similar to antibiotics.^[1]

Some necessary terminology:

The threshold concentration of bacteria that must be present in order for the virus numbers to increase is known as the proliferation threshold X_p, which is required for active therapy.
The viral inundation threshold V_I is the minimum concentration of the phage required to cause a decrease in the total concentration of viable bacteria.
The bacteria are completely cleared if the concentration of phage is greater than a certain clearance threshold V_C.

To model coevolution, some assumptions have been taken into account.

Bacterial growth follows the logistic model^[1] in the absence of phages due to the medium being a finite nutrient source.
The phage-bacteria interaction is first order with respect to both the phage and the bacteria concentrations.^[2]
There is a time delay between a bacterial cell getting infected and its subsequent lysis, called the latent period.^[3]
Bacteria can undergo mutations and become resistant to the phages. Similarly, the phages can mutate to affect the resistant bacteria. This mutation rate is proportional to the number of bacteria, as each bacteria is equally likely to mutate. The phage mutation is taken into account by assuming that the resistant bacteria mutate back to being susceptible.

The equations are as follows:

Variables

	Symbol	Details
Susceptible Bacteria	SB	Concentration of bacteria susceptible to infection by phages.
Resistant Bacteria	RB	Concentration of bacteria resistant to infection by phages.
Infected Bacteria	IB	Concentration of bacteria infected by phages.
Free Virions	V	Phage concentration in the medium.
Lysed Bacteria	LB	Concentration of dead bacteria after lysis.
Time	t	¯\_(ツ)_/¯

Constants

	Symbol	Details
Growth Rate Constant	a	Rate constant for bacterial growth.
Carrying Capacity	c	Maximum concentration of bacteria supported by the medium.
Forward Mutation Rate Constant	fr	Rate constant for mutations leading to resistance.
Backward Mutation Rate Constant	br	Rate constant for mutations leading to loss of resistance.
Interaction Constant	ir	Rate constant for phage-bacteria interaction.
Burst Size	bs	Average number of phages released per infected bacterial cell on lysis.
Decay Rate Constant	dr	Rate constant for phage decay.
Latent Period	t_L	The time delay between infection and lysis.

Dosage Determination

To determine the proliferation threshold: Consider a single infected bacterium, which upon lysis, releases bs virions. These virions either degrade, or infect more bacteria. The virions are distributed among the terms

(for viral decay)

(for virus-bacteria interaction)

R₀

R₀ > 1

To determine the inundation threshold: For the viable bacterial concentration to decrease,

The clearance threshold cannot be determined analytically from these equations. However, for any given value of bacterial concentration, the value can be determined numerically using the following equations, and taking the minimum value of the solutions:

Numerical Solutions

The problem described above is a system of five differential equations, of which two are so called delayed differential equations. They contain a term that needs to be evaluated at a time-point in the past. A python script was used to simulate the growth curves numerically, using the method of finite differences. Each of the model parameters were optimised using the Levenberg-Marquadt algorithm^[4][5] using the scipy.optimize python library. Initial values for the parameters were taken by using the data found here, and a few were chosen by hand.

The code can be found here.

Experimental Data

The protocol we executed for obtaining the experimental curves is as follows:

Put 190 μl of Lysogenic Broth (LB) and 10 μl of secondary inoculum at 0.6 OD in each well of a 96 well plate.
To each of the wells, add 2 μl of varying dilutions of phages.
Put the plate into a plate reader and take OD readings for required period of time.

The experimental data can be found here.

Getting the fit

Combinations of parameters that result in providing a fit closest to the experimental curve can most easily be found by analysing for different sets of parameters systematically.

The theoretical curve was plotted using the numerical solutions.
Then a cost function was made between the theoretical and actual curve for reasonable values of parameters.
Then scipy.optimize.least_squares was used to minimise this cost and get required parameter values.

The code for fitting can be found here.

Conclusions

Here are some of the graphs obtained after fitting our model with the experimental curves. For the bacterial concentration graphs, the experimental curve is coloured blue and the theoretical fit, orange. The viral concentrations are in the order of 4 x 10⁴ pfu.

Figure 1: SB(0) = 0.1422, V(0)=0.001

Figure 3: SB(0) = 0.1422, V(0)=0.001

Figure 2: SB(0) = 0.1408, V(0) = 0.001

Figure 4: SB(0) = 0.1386, V(0) = 0.01

Variables

	Symbol	Average Value
Growth rate Constant	a	3.76 × 10^-3
Interaction Constant	ir	6.3 × 10^-2
Burst Size	bs	1.57 × 10²
Forward Mutation Rate Constant	fr	8.1 × 10^-1
Backward Mutation Rate Constant	br	2.65 × 10^-2
Decay Rate Constant	dr	4.7 × 10^-3
Latent Period	t0	6.3 × 10²
Carrying Capacity	c	5.6 × 10^-1

Calculated values of the important dosage parameters:
Proliferation threshold X_p = 11.1
Inundation threshold V_I = 6 × 10^-2

The conditions required for active therapy could not be attained in this system, as the carrying capacity of the system was lesser than the proliferation threshold.

Limitations

The curves predicted by the model agreed with the experimental curves only upto a limit, after which there was a deviation from the ideal model. This is because of a few factors:

The model did not account for the OD of lysed bacteria, which significantly increases with time.
In the model, the backward rate of mutation is kept constant, effectively making the model similar to predator-prey dynamics. In reality, the backward rate of mutation is non zero in the beginning, but tends to 0 after some time. This accounts for the experimental curves showing no bacteria-phage interactions over large periods of time.

However, within reasonable time frames, the model predicts the phage-bacteria dynamics quite effectively, making it a useful tool to draw therapeutic correspondence between antibiotics and phage therapy.

Mutagenesis

To verify the feasibility of APES, a model was developed in parallel with the design for the hardware. The model was intended to predict the number of mutant phages produced that could infect and lyse the resistant bacteria.

To model the mutagenesis, the following assumptions were taken into account:

Each nucleotide in the genome mutates independently, with the probability of mutation following an exponential distribution.
The process is assumed to happen discretely, with mutation from one nucleotide to any nucleotide being equally likely.

The Transition Matrix

Let the rate constants for each transition i->j be given by q_ij where both i and j belong to the set {0, 1, 2, 3}. (The state 0 represents the nucleotide with base A, 1 represents base T, 2 represents base G, and 3 represents base C. The base being in the state i implies that the base is currently the base which is repesented by i) We have the constraints

Let p_ij be defined as

These definitions allow a transition matrix to be defined for the mutation, and to treat it as a discrete Markov process, with a step size given as Δt, a sufficiently small time unit. The transition matrix for the mutations, T, is now given by

where m_ij represents the element in the i^th row and j^th column of the matrix. After n steps, the probability of transition from a state i to a state j is now given by the the element in the i^th row and j^th column of Tⁿ.

The Calculations

The length of the bacteriophage T4 genome is by 168903 base pairs, out of which about 65% is assumed to be essential for the viability of the phage.^[7] Any change in any amino acids in the viable sites is assumed to result in loss in viability. The nucleotides in the genome encoding for the region of the tail fibres responsible for binding to a receptor on the bacterial cell surface are found to be in the range 156542-159622 in the T4 genome.^[8]

Given the time of mutation t, the number of steps is given by t/Δt, and hence, the transition matrix is given by Tⁿ. This is numerically computed using python and the numpy library. Using the diagonal elements of the transition matrix, which give the probability of no change in a nucleotides. These probabilities raised to the power of the number of nucleotides gives the probability that a cell remains viable. Next, the probability of transition of the nucleotide sequence between 156542-159622 to the random amino acid sequence chosen is also calculated, using data derived from an amino acid codon chart and the transition matrix. The events that there is no loss in viability and that a mutant develops which is capable of infecting the resistant bacteria are assumed to be independent. Hence, the probability of obtaining a required mutant phage is given by the product of these two probabilities. The expected value of the number of desired phages is found by multiplying the total number of phages subjected to mutagenesis times the probability calculated above.

The code used to perform the numerical analysis can be found here.

Conclusions

For a time t=500 s, and q_ij=10^-7, the probability of the virus remaining viable was found to be 10^-49. The largest number of phages used for the experiments were of the order of 10¹⁰ pfu, which gives an expected value of 10^-39 desired phages. This is an astronomically small number, clearly showing it is unreasonable to expect any kind of result from the APES experiments. Increasing the values of q_ij=10^-7, that is, using a more effective mutagen, fails to help as well, since it only results in lower amounts of viability. Decreasing the values of q_ij=10^-7 results in effectively no mutations, once again failing to produce results.

Thus, our modelling showed that no reasonable value of the mutation rates could produce desirable results. This was a major motivation for the team to decide to halt the APES experiments.

References

1]Payne, Robert JH, and Vincent AA Jansen. "Understanding bacteriophage therapy as a density-dependent kinetic process." Journal of Theoretical Biology 208.1 (2001): 37-48.
2]Jason, A. C. "A deterministic model for monophasic growth of batch cultures of bacteria." Antonie van Leeuwenhoek 49.6 (1983): 513-536.
3]Levin, Bruce R., Frank M. Stewart, and Lin Chao. "Resource-limited growth, competition, and predation: a model and experimental studies with bacteria and bacteriophage." The American Naturalist 111.977 (1977): 3-24.
4]Levenberg, Kenneth. "A method for the solution of certain non-linear problems in least squares." Quarterly of applied mathematics 2.2 (1944): 164-168.
5]Marquardt, Donald W. "An algorithm for least-squares estimation of nonlinear parameters." Journal of the society for Industrial and Applied Mathematics 11.2 (1963): 431-441.
6]Koskella, Brockhurst. "Bacteria–phage coevolution as a driver of ecological and evolutionary processes in microbial communities." US National Insitute of Medicine (2014).
7]Dedrick, Rebekah M., et al. "Functional requirements for bacteriophage growth: gene essentiality and expression in mycobacteriophage G iles." Molecular microbiology 88.3 (2013): 577-589.
8]Bartual, Sergio G., et al. "Structure of the bacteriophage T4 long tail fiber receptor-binding tip." Proceedings of the National Academy of Sciences 107.47 (2010): 20287-20292.