Difference between revisions of "Team:IISER Pune/Modeling"

Mathematical modelling suffices to – 

Provide an objective view towards complex & apparently diverse systems. This objectivity is instrumental for efforts motivated by the idea of generalization.

Provide an ability to predict the behavior of a system i.e. to predict the properties of a state of the system given the information for some other state.

Simulations –

Give us a way in which the numbers mimic biological/chemical ‘characters’ & equations derive their behavior.

Help us in visualizing the results of equations that are analytically unsolvable.


Variable link oscillator

Description of original model from reference[1]

This is a quantitative model describing the regulation of the PRM operator region of λ phage. The model envisions that the system is a DNA plasmid consisting of the promoter region and cI gene. 
The promoter region contains the three operator sites known as OR1, OR2, and OR3. 
The gene cI expresses repressor (CI), which in turn dimerizes and binds to the DNA as a transcription factor. This binding can take place at one of the three binding sites OR1, OR2, or OR3. The binding affinities are such that, typically, binding proceeds sequentially; the dimer first binds to the OR1 site, then OR2, and lastly OR3. Positive feedback arises due to the fact that downstream transcription is enhanced by binding at OR2, while binding at OR3 represses transcription, effectively turning off production and thereby constituting a negative feedback loop. 
The chemical reactions describing the network are naturally divided into two categories – fast and slow. The fast reactions have rate constants of order seconds, and are therefore assumed to be in equilibrium with respect to the slow reactions, which are described by rates of order minutes. If we let X, X2, and D denote the repressor, repressor dimer, and DNA promoter site, respectively, then we may write the equilibrium reactions
  

where Di denotes dimer binding to the ORi site, and the Ki = ki/k−i are equilibrium constants. We let K3 = σ1K2 and K4 = σ2K2, so that σ1 and σ2 represent binding strengths relative to the dimer-OR1 strength. 
The slow irreversible reactions are transcription and degradation. If no repressor is bound to the operator region, or if a single repressor dimer is bound to OR1, transcription proceeds at a normal unenhanced rate. If, however, a repressor dimer is bound to OR2, the binding affinity of RNA polymerase to the promoter region is enhanced, leading to an amplification of transcription. Degradation is essentially due to cell growth. We write the reactions governing these processes as 
 
where P denotes the concentration of RNA polymerase, n is the number of repressor proteins per mRNA transcript, and α > 1 is the degree to which transcription is enhanced by dimer occupation of OR2. 
Defining concentrations as our dynamical variables, x = [X], x2 = [X2], d0 = [D], d1 = [D1], d2 = [D2D1], and d3 = [D3D2D1], one can write a rate equation describing the evolution of the concentration of repressor, 
 
where they assume that the concentration of RNA polymerase p0 remains constant during time. 
They next eliminate x2 and the di . The model utilizes the fact that the reactions in Eq. (1) are fast compared to expression and degradation, and write algebraic expressions 
Further, the total concentration of DNA promoter sites dT is constant, so that

 
where 
m is the copy number for the plasmid, i.e., the number of plasmids per cell.
They eliminate two of the parameters by rescaling the repressor concentration x and time. 
Now, they define the dimensionless variables x􏰄’ = x (K1K2)1/2  and t’ = t(ktp0dT n(K1K2)1/2). Upon substitution, 
 

The parameter γx is directly proportional to the protein degradation rate, and in the construction of artificial networks, it can be utilized as a tunable parameter. The integer parameter m represents the number of plasmids per cell. While this parameter is not accessible during an experiment, it is possible to design a plasmid with a given copy number, with typical values in the range of 1-100. 
The hysteretic effect can be employed to induce oscillations, provided we can couple the network to a slow subsystem that effectively drives the parameter γx. This can be done by inserting a repressor protease under the control of a separate PRM promoter region. On one plasmid, we have the network of the repressor protein CI, which is under the control of the promoter PRM , stimulates its own production at low concentrations and shuts off the promoter at high concentrations. On a second plasmid, again utilize the PRM promoter region, but here they insert the gene encoding the protein RcsA. The crucial interaction is between RcsA and CI; RcsA is a protease for repressor, effectively inactivating its ability to control the PRM promoter region.
First, both RcsA and repressor are under the control of the same promoter, so that the functional form of the production term f(x)  will be the same for both proteins. Second, they envision the network as being constructed from two plasmids – one for repressor and one for RcsA, and that we have control over the number of plasmids per cell (copy number) of each type. Lastly, the interaction of the RcsA and repressor proteins leads to the degradation of repressor. Putting these facts together, and letting y denote the concentration of RcsA, we have 
  
 We simulated this model, in python using following set of parametric values – 
 
to get –
 


 
[C1] = 7.7x + 3.0x2

We then ‘played around’ with initial conditions –

 


Inference – Initial conditions do not affect the oscillatory period, when plasmid copy number ratio is kept constant.

  The next step was to check the effect of mx : my ratio on the oscillatory system –
 
Here the ratio is kept constant – 10:1 throughout, while the values are being changed. 
Inference -  No change in the oscillatory period. However, the concentration amplitudes at which the system oscillates exhibit differences.
We then simulated the model for a range of values for mx : my ratio – 

 
Here the initial conditions are constant for all values of mx : my ratio.
Inference – For the given parametric values (referred from the paper),  the oscillatory regime is possible only for a range of plasmid copy number ratios sarting at 4:1 and giving best results (oscillations of period reuired by us) at close to 10:1.
According to our project idea, it is crucial for the entire device to be packaged in one plasmid. However, this would mean –
   
mx = my
We came up with a no. of possible mitigation strategies / alternetives for this –
1.Multiple promoter sites (i.e. gene sets)-
	Concerns : 
	Not “effective” as plasmid copy no. ratio
	Issue of promoter titration
	Larger plasmids required
2. Different RBS Efficiency ratios-
	Advantages :
	Present as bio-bricks (measurements performed)
	Efficiency ratios comparable in magnitude to copy no. ratio (from 4:1 to         10:1)
	
	
Bio-brick part	Efficiency
BBa_B0034	1.0 (Reference)
BBa_B0030	0.6
BBa_B0032	0.3
BBa_B0031	0.07
BBa_B0033	0.01

We now simulated this ‘situation’, wherein the plasmid copy number is the same for genes, but the RBS for them have different efficiencies. These RBS efficiencies are considered, if not equal, proportional to translational efficiency.
We scale  (multiply equations with appropriate factors) the equations of our oscillatory system such that –
Ratio of RBS efficiencies = Plasmid copy number ratio for which oscillations are observed.
 
RBS Efficiency Ratio 
0.6 : 0.07 = 8.75

Further work –
1.	To construct a model that includes translational & transcriptional considerations explicitly.
2.	Attempt a combinatorial way - 
	(no. of gene sets) * (RBS efficiency ratio) ≅ plasmid copy no.
	Verify the linearity!
	
	 

Terminator (Killer Gene) Model –
Description of model (from the reference)
The paper describes a circuit that programmes a bacterial population to maintain a cell density that is lower than the limits imposed by the environment (for example, by nutrient supply). This is achieved by a quorum sensing module attached to a killer gene. 
The LuxI protein of the well-characterized LuxI/LuxR system from the marine bacterium Vibrio fischeri synthesizes a small, diffusible acyl-homoserine lactone (AHL) signalling molecule. The AHL accumulates in the experimental medium and inside the cells as the cell density increases. At sufficiently high concentrations, it binds and activates the LuxR transcriptional regulator, which in turn induces the expression of a killer gene (E) under the control of a luxI promoter (pluxI). Sufficiently high levels of the killer protein (CcdB) cause cell death. The CcdB portion retains the toxicity of native CcdB, which kills susceptible cells by poisoning the DNA gyrase complex. 
A simple mathematical model is described in the reference that predicts that the system would reach a stable cell density for all realistic parameter values, although it might go through damped oscillations while approaching the steady state. 
 
They model the major kinetic events dictating the circuit function, including cell growth and death. Assumptions of the model –
 (1) without the circuit, changes in viable cell density follow logistic kinetics with an intrinsic per capita growth rate of k and a carrying capacity of Nm 
(2) for circuit-regulated growth, the cell death rate is proportional to the intracellular concentration of the killer protein (E) with a rate constant of d (h)
(3) the production rate of E is proportional to AHL concentration (A, assumed to be the same inside and outside the cells) with a rate constant of kE
(4) AHL production rate is proportional to N with a rate constant of nA 
(5) degradation of the killer protein and AHL follows first-order kinetics with rate constants of dE and dA 

  Analytical solution for this model :
 
Our work :
We simulated this model in Python, for given values of parameters & constants –
(These values were referred from the paper)
 











 

Here, 

NC = population density when QST circuit is ‘OFF’
 N = population density when QST circuit is ‘ON’

 




We also checked the behavior of A & E with respect to N.

We then came up with an idea, where we could delay the onset of terminator.
This is crucial for our device idea, as the terminator should kick in only after the system has been hijacked to extent that it produces pigment enough to be recognized visibly.
# add schematic

Hence, we coupled our detection module with the terminator module. Depending on the whether the concentration of pigment has crossed a threshold value, the killer gene is activated.

Pigment production gene under a constitutive promoter, accumulates continuously. When sufficient (visible to eye, here assumed to be 1/10 conc. of that of lycopene in tomato) amount of pigment is accumulated, the sender (AHL producing part of the circuit) of the terminator is induced, which then in turn induces the killer gene in the receiver.
This pigment- sender induction can either be manual (as this model now helps us to know the time period in which ‘sufficient’ amount of pigment accumulates) or through a density counter associated with the pigment production gene.
   

 

N= QST on
Nc=QST off
Nd=Delayed initiation

References :
1)	Designer gene networks towards funamnetal cellular control, Hasty et al, Chaos (2002)
2)	Programmed population control by cell-cell communication and regulated killing, You et al, Nature (2004)