In 1967 Fredrickson et al. studied mathematically development of a bacterial population, under the assumptions of a large population of independant bacteria in a well mixed solution of constant volume. The large population ensures that for the population the expectation value is a good estimate of the average. The bacteria being independant ensures that the behaviour of each individual depends only on its internal state <math>\mathbf{z}</math> and the conditions <math>\mathbf{c}</math> which are the same for all individuals. The volume is well mixed so the conditions <math>\mathbf{c}</math> which are the same everywhere. The volume is constant so that the population caracteristics can be evaulated by integration over the volume. In their development <math>\mathbf{z}</math> and <math>\mathbf{c}</math> are considered to be arbitrary vector quantities.
From this starting point they develop a pair of master equations of change to describe the evolution of the population: <math>\begin{gathered}
\frac{\partial}{\partial t} W_{\mathbf{Z}}(\mathbf{z},t) + \nabla_{\mathbf{Z}} \cdotp [(\mathbf{\beta} \cdotp \bar{\mathbf{R}}(\mathbf{z,c}) W_{\mathbf{Z}}(\mathbf{z},t)] \\ = 2 \int \sigma (\mathbf{z',c}) p(\mathbf{z,z',c}) W_{\mathbf{Z}}(\mathbf{z'},t)\mathrm{d}v' - (D+\sigma (\mathbf{z,c})) W_{\mathbf{Z}}(\mathbf{z},t)\end{gathered}</math> <math>\frac{d\mathbf{c}}{dt} = D(\mathbf{c_f} - \mathbf{c} ) + \mathbf{\gamma} \cdotp \int \bar{\mathbf{R}}(\mathbf{z,c}) W_{\mathbf{Z}}(\mathbf{z},t)\mathrm{d}v</math> In these equations the various symbols are as follows:
<math>\mathbf{z} </math> | Vector for internal state of a bacteria. |
<math>\mathbf{c} </math> | Time dependant vector for conditions. |
<math>W_{\mathbf{Z}}(\mathbf{z},t) </math> | Distribution of bacteria in <math>\mathbf{z}</math> space at time <math>t</math>. |
<math>\bar{\mathbf{R}}(\mathbf{z,c}) </math> | The expected value or the reaction rate vector function of <math>\mathbf{z}</math> and <math>\mathbf{c}</math>. |
<math>\sigma (\mathbf{z',c}) </math> | Rate of fision for bacteria as a scalar function of <math> \mathbf{z,c} </math>. |
<math>p(\mathbf{z,z',c}) </math> | Partitioning probability of generating a child in state <math>\mathbf{z}</math> from a parent in state <math>\mathbf{z'}</math> given the conditions <math>\mathbf{c}</math>. |
<math>\nabla_{\mathbf{Z}}\cdot \mathbf{V} </math> | <math>\sum \frac{\partial}{\partial z_i}\mathbf{V}_i </math> |
<math>\textrm{d}v' </math> | Integral over state space <math>v'</math> . |
<math>D </math> | Dilution rate of the culture (for fermenters). |
<math>\beta </math> | Stochiometric matrix for cellular substances. |
<math>\gamma </math> | Stochiometric matrix for extra-cellular substances. |
With these relations:
<math>\bar{\dot{\mathbf{V}}}(\mathbf{z,c})</math> | <math>= \mathbf{\beta} \cdotp \bar{\mathbf{R}}(\mathbf{z,c}) </math> | The expected internal state change rate vector. |
<math> -\mathbf{\gamma} \cdotp \bar{\mathbf{R}}(\mathbf{z,c}) </math> | The expected consumation of substances in the environment by a cell in state <math>\mathbf{z}</math> . |
Thus for a particular problem in hand it is necessary to chose <math>\mathbf{z}</math> and <math>\mathbf{c}</math> that represent the state of cells and the media. Then the matrices and functions <math>\beta</math>, <math>\gamma</math>, <math>\bar{\mathbf{R}}(\mathbf{z,c})</math>, <math>\sigma (\mathbf{z',c})</math> and <math>p(\mathbf{z,z',c}) </math> need to be defined for the problem considered. Finally the inital conditions <math>W_{\mathbf{Z}}(\mathbf{z},t)</math> and <math>\mathbf{c}_0 </math> and growth conditions <math>D</math> and <math>\mathbf{c}_f </math> need to be fixed.
For the problem in hand, plasmid maintenance during growth with 2 different plasmids, and attempting to find a simple solution to the problem we propose a 3 variable internal state vector:
<math>\mathbf{z} = \begin{bmatrix} z_0 \\ z_1 \\ z_2 \end{bmatrix} =
\begin{bmatrix}\textrm{Cell maturity} \\ \textrm{count of plasmid 1} \\ \textrm{count of plasmid 2} \end{bmatrix}</math>
In this internal state vector: <math>z_0</math> is a mesure of the growth of the bacteria, encompassing such things as size, number of chromosomes and mass; <math>z_1</math> and <math>z_2</math> represent the number of copies of each plasmid. For the external conditions we propose simply the substrate concentration <math>S</math>. The maturity parameter has a minimum value of 1 and must increase to 2 before division can occur.
For the rates of change of the internal state vector then we propose for the bacterial maturity to extend the development presented in Shene et al. to include 2 plasmids and incorporate the cell maturity as a state variable. This gives: <math>\dot{z}_0 (\mathbf{z},S) = \mu = \mu _{max} \frac{S}{K_S+S} \frac{K_{z_1}}{K_{z_1}+z_1^{m_1}} \frac{K_{z_2}}{K_{z_2}+z_2^{m_2}}</math> Here <math> \mu _{max}</math> is the maximum growth rate <math>hr^{-1}</math>: <math>\mu (\mathbf{z,S})</math> the growth rate ; <math>K_S</math> is the Monod constant in <math>g/l</math> for the substrate; <math>K_{z_1}</math> is the inhibition constant for plamid number 1 in (plasmids per cell)<math>^{m_1}</math>, and <math>m_1</math> the Hill coefficient for the cooperativity of inhibition. <math>K_{z_2}</math> and <math>m_2</math> represent the same parameters for plasmid 2.
For plasmid replication rate we propose, again following Shene et al. , the empirical relationship : <math>\dot{z}_1 (\mathbf{z},S) =
\begin{cases} k_1 \frac{\mu (\mathbf{z},S)}{K_1 + \mu (\mathbf{z},S)} ( z_{1_{max}} - z_1 ) & \text{if } z_1 \geq 1.0 \\ 0 & \text{if } 0.0 \leq z_1 < 1.0 \\ \end{cases}</math> This relation, and equivalent one for plasmid number 2 <math>\dot{z}_2 (\mathbf{z,S})</math> is designed to satisfy the boundary conditions of no reproduction if there is less than 1 plasmid in the cell, and a maximum copy number of <math>z_{1_{max}}</math>. This introduces the parameters <math>k_1</math> and <math>K_1</math> which are respectively the plasmid replication rate (in <math>hr^{-1}</math>) and the inhibition constant (also in <math>hr^{-1}</math>). The inhibition constant reduces plasmid replication rate at slower growth rates.
Notice that here we have directly defined <math>\bar{\dot{\mathbf{V}}}(\mathbf{z,c})</math> rather than <math>\beta </math> and <math>\bar{\mathbf{R}}(\mathbf{z,c})</math>. For the growth yield we propose : <math>\gamma \cdotp \bar{\mathbf{R}}(\mathbf{z,c}) = \alpha \mu(\mathbf{z},S)</math> Where alpha is the growth yield in <math>g/l/cell</math>. The remaining functions and parameters in equations 1 and 2 are the division rate <math>\sigma (\mathbf{z,c})</math> and the partitioning function <math>p(\mathbf{z,z',c}) </math>. There is less consensus in the litterature for an at least empirically appropriate form for these equations. To remain simple we propose: <math>\sigma (\mathbf{z,c}) = \sigma \times H[2.0] =
\begin{cases} 0 & \text{if } z_0 < 2.0 \\ \sigma & \text{if } z_0 \geq 2.0 \end{cases}</math> Here we assume that there is a fixed rate of division <math>\sigma </math> once cells are big enough to divide (<math>H[]</math> is the Heaviside function). <math>p(\mathbf{z,z',c}) = p(z_0,z'_0) \times p(z_1,z'_1) \times p(z_2,z'_2)</math> <math>p(z_0,z'_0) = \delta_{z_0,\frac{z'_0}{2}} = \begin{cases} 1 & \text{if } z_0 = z'_0/2.0 \\ 0 & \text{if } z_0 \ne z'_0/2.0. \end{cases}</math> <math>p(z_1,z'_1) = 0.5^{z'_1} \begin{pmatrix} z_1 \\ z'_1 \\ \end{pmatrix} = 0.5^{z'_1} \frac{z'_1 !}{(z'_1-z_1)! z_1 !}</math> In these equations we assume that the partitioning of the three internal state variables are independant. That cells divide exactly in half, that is the maturity parameter is exactly halved when the cells divide (<math>\delta</math> is a Kronecker delta function). That the two plasmids segregate independantly and as individual plasmids according to a binomial distribution. These assumptions are probably the most suspect in the model.
This initial version of the model has no contention, that is <math>z_1</math> and <math>z_2</math> do not influence the growth rate <math>\mu </math>. In order to develop the model for the system envisaged this needs to be introduced.
cccp0.55 Symbol & Value & Units & Comments
<math>S_f</math> &10.0& <math>g/l</math> & Substrate concentration in feed.
<math>S_0</math> &10.0& <math>g/l</math> & Initial substrate concentration in fermenter.
<math>D</math> &0.5& <math>hr^{-1} </math>& Dilution rate
<math>\alpha </math> &0.5& <math>g/l/cell</math> & Growth yield on substrate.
<math>\mu _{max} </math> &2.0& <math>hr^{-1}</math> & Maximum growth rate on substrate.
<math>K_S</math> &0.25& <math>g/l</math> & Monod constant for substrate.
<math>K_{z_1}</math> &100.0& <math>(cell^{-m_1})</math> & Plasmid 1 growth inhibition constant
<math>m_1</math> &1& & Plasmid 1 Hill coefficient
<math>k_1</math> &20& <math>hr^{-1}</math> & Plasmid 1 replication rate.
<math>K_1</math> &0.1& <math>hr^{-1}</math> & Plasmid 1 replication inhibition constant.
<math>z_{1_{max}} </math>&10&& Plasmid 1 maximum copy number.
<math>K_{z_2}</math> &100.0& <math>(cell^{-m_1})</math> & Plasmid 2 growth inhibition constant
<math>m_2</math> &1& & Plasmid 2 Hill coefficient
<math>k_2</math> &20& <math>hr^{-1}</math> & Plasmid 2 replication rate.
<math>K_2</math> &0.1& <math>hr^{-1}</math> & Plasmid 2 replication inhibition constant.
<math>z_{2_{max}} </math>&10&& Plasmid 2 maximum copy number.
<math>\sigma</math> &5.0& <math>hr^{-1}</math> & Rate constant for cell division once cells are big.
<math>k_{a1} </math> &2.0&& Toxin 1 efficiency parameter
<math>k_{a2} </math> &0.5&& Toxin 2 efficiency parameter
<math>k_b</math> &2.0&& Ratio of anti-toxin to toxin production rates
Substituting into the equations 1 and 2 we obtain: <math>\begin{gathered}
\frac{\partial}{\partial t} W_{\mathbf{Z}}(\mathbf{z},t) + [\nabla_{\mathbf{Z}} \cdotp \bar{\dot{\mathbf{Z}}}(\mathbf{z},S)] W_{\mathbf{Z}}(\mathbf{z},t) + \sum_i \bar{\dot{z_i}} \times \frac{\partial}{\partial z_i} W_{\mathbf{Z}}(\mathbf{z},t) \\ = 2 \sigma \int_{z'_0>2.0} \delta_{z_0,\frac{z'_0}{2}} \times 0.5^{z'_1+z'_2} \times \begin{pmatrix} z_1 \\ z'_1 \\ \end{pmatrix} \times \begin{pmatrix} z_2 \\ z'_2 \\ \end{pmatrix} \times W_{\mathbf{Z}}(\mathbf{z'},t)\mathrm{d}v' \\ - (D+\sigma H[2.0] W_{\mathbf{Z}}(\mathbf{z},t))\end{gathered}</math> <math>\frac{d\mathbf{c}}{dt} = D(c_f - c ) - \alpha \int \mu (\mathbf{z},S) W_{\mathbf{Z}}(\mathbf{z},t)\mathrm{d}v</math>
The aim of studying the behaviour of this model is to investigate how growth conditions <math>D, S_f</math> and time modulate the development of the population in state space <math>W_{\mathbf{Z}}(\mathbf{z},t)</math>. In particular we are interested in finding how the number of bacteria without plasmids <math>W_{\mathbf{Z}}([z_0,0,0],t)</math> in the culture progresses, and how this depends on the various parameters in the model.
We need to introduce a modification to equation 3 (the state dependant growth rate) in which the equilibrium between <math>z_1</math> and <math>z_2</math> features. This modification is to take into account that each toxin should be inhibited by anti-toxin for maximal growth. A possibility is that we consider that for each toxin anti-toxin pair:
This model of toxin anti-toxin interaction takes into account the bacteriostatic nature of most such toxins, and the presence of measurable <math>IC_{50}</math> values. Clearly a more realsitic model would need to take into account cell volume, protein synthesis rates etc. Nevertheless, this simple model gives:
<math>\begin{aligned}
\mu &= \mu_{0} \times e^{-\frac{[Toxin]}{IC_{50}}} \\ [Toxin] &= max(0,\frac{k_1z_1-k_2z_2}{k_3}) \\ \mu &= \mu_{0} \times min(1.0,e^{-k_a(z_1-k_bz_2)}) \\ k_a &= \frac{k_1}{k_3 \times IC_{50}} \\ k_b &= \frac{k_2}{k_1} \end{aligned}</math> Incorporating 2 toxin anti toxin pairs, if we assume that the production rate ratio is the same for both, <math>k_b</math>, is independant of the system we have: <math>\mu = \mu_{0} \times min(1,e^{-k_a(z_1-k_bz_2)})) \times min(1,e^{-k_a(z_2-k_bz_1)}))</math> where <math>\mu_0</math> is given by equation 3. For each toxin the efficiency parameter is a measure of ratio of toxin accumulation in cells with one gene copy and without anti-toxin to the <math>IC_{50}</math>.
The equations are hard to integrate and with partial derivatives did not seem easy to tackle with standard packages. We therefor wrote a simple custom integrator that calculates the growth trajectory of about <math>10^6</math> bacteria in a fermenter at constant dilution rate. This essentially implements a Monte-carlo Markov-chain integrator.
The first version of this program is written in Python, and is shown and described here.
The first lines of the program load various modules that are used, these are: “math”, required for the <math>exp()</math> function; “copy”, for duplicating a list of individuals with the <math>deepcopy()</math> function, and the numerical package “numpy”, for the random number generator.
A simple function “michaelis” is then defined. As several of the model equations contain elements in the form of the Michaelis Menten equation, for example equations 3 and 4, this makes the code easier to read.
The program then initializes the various model parameters. These include the parameters in table [tbl:parameters], and a few other parameters specific to the integrator. These extra parameters are: <math>V</math>, the observation volume in ml; <math>V_{tot}</math>, the total fermenter volume (taken as 100 ml); <math>total</math>, the total number of bacteria in the integration volume <math>V</math>; <math>t_{max}</math>, the total time to run the integrator for, in hours; <math>dt</math>, the time step, in hours. In addition the dilution parameter, <math>D</math>, is changed from a simple first order rate constant to be a flow rate in ml/hr. This can be divided by <math>V_{tot}</math> to obtain the first order rate constant. In this section the time counter and the random number generator are both initialized.
An empty list of bacteria, “population” is created. Various integration parameters are zeroed. The next lines create the initial population of cells, assigning a random growth state between 1 and 2 to each cell, and giving the maximum number of plasmids. It was decided to give each cell the maximum number of plasmids, this was done to avoid the rapid creation of a significant population of plasmid free cells, due to the division of cells with few plasmids early in the simulation, when this quantity was randomized.
In the next lines of the program we enter the integration loop which runs until the timer expires. The first operation is output the current state of the population, line 92. Then the various counters that keep track of the population status are re-initialized ready for the next time-step. At line 103 the program enters the loop to update each bacterium in the population. This loop runs backwards (from the end to the begining of the list), so we can delete the current bacterium if necessary. we then collect a random number to determine the fate of the current bacterium. If this number is less than the probability of a bacterium being lost from the simulated volume <math>V</math> in a timestep by washout then the current bacterium is destroyed. The random value is reduced by the probability of washout so it can be used for further testing.
If the bacterium survives the washout then it grows. Here the values of <math>\dot z_0</math>, <math>\dot z_1</math> and <math>\dot z_2</math> are calculated using equations 3, 4 and 17. Then the state vector of the current bacterium is updated. In line 135 the total growth of bacteria during the time-step is summed, this is used to calculate the substrate use later, line 147.
After the cells, that avoided being washed out have grown, those that are large enough to divide (<math>z_0 > 2.0</math>) are given a chance to divide, depending on the random control parameter “cntrl”. If the cell is selected for division a new bacterium is created and allocated about half of the resources. The other half of the resources remain in the original bacterium. The allocation of the plasmids (<math>z_1</math> and <math>z_2</math>) follow a binomial distribution, equation 9. The division of the remaining cell contents, was adjusted to follow a normal distribution, with mean 0.5 the parents content and standard deviation 0.05. This was implemented rather than the simple division proposed in equation 8 and 10 to avoid persistent synchrony and oscillations in cell number during slow growth.
Having updated all the cells, including additions and subtractions due to cell division and washout, the new total number of cells is obtained and the substrate concentration updated according to equation 11.
In the last lines of the program the total cell number is controlled to avoid it becoming either too large or too small, this is achieved my modifying the simulated volume.
This Python program handles the integration as the growth of a population in fermenter like conditions, essentially following the equations set out above.
Equation 17 describes how growth rate depends on the number of each type of plasmid, and 5 of the parameters in table [tbl:parameters] namely <math>K_{z1}</math>, <math>K_{z2}</math>, <math>k_{a1}</math>, <math>k_{a2}</math>, and <math>k_{b}</math>. The graph of this function growth rate (<math>\mu </math> or <math>\dot Z_0</math>), as a function of the number of each type of plasmid (<math>Z_1</math> and <math>Z_2</math>), is shown in figure [fig:contour] as a contour plot. This shows a central flat surface sloping away from the no plasmid situation, that is cut away close to the axes where the amount of anti-toxin is no longer enough to prevent toxicity. In the parameter space the slope of the central surface is determined by <math>K_{z1}</math> and <math>K_{z2}</math>, the more expensive plasmid maintenance is the faster the fall off. The width of the central wedge is determined by the parameter <math>k_b</math> larger values give a wider wedge, increased anti-toxin production reduces toxicity due to plasmid imbalance. The fall off in growth rate outside this central region depends on the toxicity parameters <math>k_{a1}</math> and <math>k_{a2}</math>.
We can this note that the best growth rate with plasmids is with one of each, and no plasmids is always best. This is because in our model <math>K_{z1}</math> and <math>K_{z2}</math> are positive. As mutations, which are not included in the model, will tend to increase growth rate they will tend to increase <math>k_b</math> and decrease the other parameters. To minimize the risk of plasmid loss the parameter <math>k_b</math> would seem particularly critical. If the wedge is too narrow all plasmids become toxic, if it is too wide growth with few of one plasmid becomes possible because of insufficient toxicity.
As the model for plasmid segregation is a binomial distribution, with probability of 1/2 for each daughter cell, the probability of generating a daughter with no plasmids in a division is <math>2 \times \frac{1}{2^n}</math>, where n is the total number of plasmids to segregate. The prefactor of 2 is as each division gives 2 daughter cells. Thus with 2 plasmids with a copy number of 50 the probability of generating a plasmid free daughter is <math>2^{-99}</math> or about <math>10^{-30}</math>. As a 1.0 <math>OD_{600}</math> culture corresponds to about <math>10^8</math> cells/ml in 100ml of culture growing at constant volume at OD=1.0 and a doubling rate of 2 <math>hr^{-1}</math> the probability of this arriving is about 1 in <math>10^{20}</math> per division, or once in <math>5 \times 10^{19}</math> hours.
However the rate is likely to be increased as not all cells will have 100 plasmids when they divide, and the probability is very sensitive to this. With only 20 of each plasmid/cell the frequency shoots up, under the same growth conditions, to 1% of the times the <math>10^{10}</math> cells divide, ie on average once in a hundred generations or every 2 days.
Therefor critical to the maintenance of the plasmids is the number of plasmid copies in dividing cells. This needs to be evaluated by integration.
0.8CCCC D(ml/hr)&S(g/l)&Cell density(<math>OD_{600}</math>)&Average number of plasmids
5 & 0.0030 & 2.05 & 19.87
10 & 0.0063 & 2.02 & 19.79
20 & 0.0138 & 1.99 & 19.60
50 & 0.0473 & 1.89 & 18.99
100 & & &
200 & & &
250 & 10.00 & 0.00 & n.a.
Simulations were run using the model to see how the cell density and substrate concentration adjusted to the dilution rate. As the number of plasmids in cells is an important parameter this was also calculated. These results, shown in table [tbl:density], demonstrate how the equilibrium cell density decreases as dilution rate increases, while the equilibrium substrate concentration increases. The average number of plasmids per cell also decreases as growth rate increases. In these simulations, as the maximal number of each type of plasmid was 10, see table [tbl:parameters], plasmid loss was a frequent occurence. As can be seen once the flow rate gets too great the above about 250 ml/hr the bacteria are washed out and a stable culture is not formed.
TODO: next experiments number of plasmids at division for 100 and 200 with max = 20 (graph) p(n) vs n for 2 conditions.
TODO: then calculate prob child with 0 as <math>\int P(n) 0.5^(n-1) dn</math> for a dividing bacterium.