Team:HiZJU-China/Model

• Home
  Homepage
• Project
  Description Design Result Hardware Proposed Implementation Proof of Concept
• Team
  Members Attributions
• Lab
  Experiments Notebook Protocols Model Safety
• Parts
  Parts Engineering Success Contribution
• Human Practices
  Overview Intergrated Human Practice Education and Public Engagement Collaborations
• Awards
  Judging Sustainable Development Impact Best Safety

PartⅠ-Degradation Modules

Model Introduction
TTo trace the variables in the whole process of detection and decomposition, we decide to build a mathematic model to give out an illustration. By using some simple theories, the model was established and well presented.

Goals
1. Explaining the construction of the detection model and using this to prove that the detection approach is still work under a competitive situation.
2. Illustrating the relationship between the variables and the time in the decomposition. Analyzing the sensitivity of parameters in the reaction and the efficiency of reaction.

Theories
Basically, we use several kinetic laws to construct the model.
1. Mass action
In chemistry, mass action is that the rate of a reaction is proportional to the concentration of reactants. For example, the forward rate of following equilibrium

is

And that will lead to the result that when this system is stable, the proportion of each substance will be a constant, which can be calculate as follow:

2. Hill equation
To illustrate the proportion of the ligand being combined to the recipient, Hill provided a method to quantify. It can be presented as follow. Suppose we know the K and the total concentration of ligand. Then the fraction is

However, this uses an assumption that the recipient is considered as constant.

3. Michaelis-Menten equation
Michaelis and Menten’s study points out that an enzyme reaction can be seen as this:

And when the reaction is in a stable state, the concentration of ES shall be a constant. Then we can naturally get the rate of the starting point is

where Ks is Kr/Kf .
However, this is imprecise because the second step of the reaction can still be reversible. Following this, we immediately get the advanced version of the equation.

And the rate is

where Km is (Kr1+Kf2)/Kf1 .
If we suppose that the reversible equation is right, then the rate during the reaction can be presented as

where k1 is kf1, k-1 is kr1 and so on so forth. This is also an approximation.

4. The ODE system
The actual rate of a substance is

If we use this to every substance in the system, we can get an ODE system to describe the concentration of any substance at any time. This is the way we get the precise answer of the quantities in the system. But comparing to the methods above, it requires more time to calculate and harder to be implemented.

The Detection
Using the Hill equation, the maximum bounded fraction is

Here [Ks]^n=K. By using this, we can predict the limitation of product.

Fig 1. Theoretic bounded proportion. Ks=kb=1.

The limitation is the proportion of bounded substrate.
However, this isn’t precise. In reality, n is not exactly equal to the number of places a ligand can bound. Hemoglobin has four places for the ligand, but its n is 3.0. This is because the ligand bounded will prevent the bound of the other ligand, and a recipient probably only has few ligands bounded, which is not full. We have to hence change the model, simply by introducing the different levels of Ks. After we have changed the model, it can be seen the difference between the two curves.

Fig 2. Changed model.

In practice, the bounded proportion of the ligand positively correlates with the amount of the fluorescent protein. But in our real detection, we have to face a problem that two ligands are competing for the same recipient. To show this, we use the simbiology in MATLAB to get a fast simulation.

Fig 3. The diagram of two ligands competing for single recipient.

Fig 4. L1 and L2 compete for the recipient.

Now the k_f1=10, k_f2=1, L1=1, L2=10, R=1, the other constant is one. Through the K1 is ten times of K2, it can’t compensate the loss of L1. If we try to let k_f1=1, k_f2=1, k_r1=0.1, where the K1 remains unchanged, it will be like this.

Fig 5. Changing part of the constants.

The final proportion of RL1 is still the same, but it takes a lot longer to enter the stable state.
Then we change k_f1=200 and get new result.

Fig 6. Increased the reaction speed.

This time RL1 is far more than RL2, and the selectivity of the recipient we want finally comes out.
So, to make the detection precise and fast, the thing we need to do is to make k_f1 as big as possible.

The decomposition
Firstly, we would like to see the theoretic limitation of a simple enzyme reaction. The reaction can be presented as

The relationship between the concentration of product and the time can be figured out by using ODE system. If we use the reversible Michaelis-Menten equation to approximate the rate during the reaction, the differential equation is

Solving this equation, we can get

While we can’t get the exact form of P(t), we do get the inversed version, which is t(P) presented above. This lead to the first graph.

Fig 7. Reversible Michaelis-Menten equation.

However, this is not ideal for the concentration of ES is not guaranteed to be constant throughout the reaction. And also the substrate can’t be reacted thoroughly. Therefore, we develop a bigger differential equation system to solve this. To take a smaller step, we first assume that the substate is constant. Using the theory of mass action, the equations are follow.

The curve can then be compute.

Fig 8. Approximate ODE.

This looks similar to the curve in the figure 3 and reveals the partial correctness of using Michaelis-Menten equation to predict the whole procedure. To get the final answer, we discard the promise that the substrate is constant. This result to add only one new equation to the system, and it finally looks like below.

Fig 9. Precise ODE.

Putting them altogether and use hill equation to prove the correctness, we can get

Fig 10. Altogether.

Here k_1=1, k_(-1)=10, k_2=1, k_(-2)=10. See that only the final system can illustrate the reaction corresponding to hill’s prediction.

And now we can finally build our decomposing system. This involves 4 reactions.

Then we can use simbiology to illustrate it.

Fig 11. Figure of the decomposition system.

The system is competitive inhibition, which means NO and OH-EE2 are both competing for the HAO and AMO at the same time. So enhancing the rate of the two enzyme reactions can’t necessarily increase the proportion of the OH-EE2 in the product.

Fig 12. Kf1 equals 20.

Fig 13. Kf1 equals 0.1.

As you can see, the only thing the change of Kf1 affects is the total reaction rate. The concentration of the product depends on the non-enzyme reactions. We have done the sensitivity calculation to prove this.

Fig 14. Sensitivity calculation.

However, these non-enzyme reactions can’t be changed manually, so we can only try to assume the states during the reaction. Four possibilities may happen, here we only present two of them.

Fig 15. OH-EE2 is produced fast but ends up less.

Fig 16. OH-EE2 is produced slow and ends up more.

To figure out which one is correct, we have done several experiments. Taking the real data into our model, this is what we get.

Fig 17. Approximation of the real data. Initial concentration of EE2 is 5.06×10^(-5) mol/L.

Because of the limitation of the computer and the experiment data, we can give out the whole precise graph. But the information from this graph can be translated as after 6 hours of decomposition, the concentration of the EE2 has dropped about one eighth. This efficiency is acceptable. Further exploration may need to get a better result.

PartⅡ-Detection Modules

Abstract: In order to reduce the detection limit of the biological detection system, we modified the ligand binding domain (LBD) of estrogen receptor by LeDock prediction to improve the specificity and affinity of receptor. Then PLIP was used to get a better understanding of mutations. The following processes and methods led us to explore the optimal solution.

Introduction to estrogen receptor ERα
Estrogen receptors are members of the superfamily of nuclear receptors, which are responsible for mediating the effects of estrogens in target tissues by acting as ligand dependent transcription factor. ER is composed by three structural domains: a modulating domain with ligand-independent transactivation, a DNA-binding domain (DBD) and a ligand binding domain (LBD).

There are two subtypes of human estrogen receptors (ERα and ERβ) that differ in tissue distribution. It has been reported that when transcription from a gene depends on both AF1 and AF2, the activity of ERα greatly exceeds that of ERβ. To order to enhance the detection ability of EE2, we choose ERα as biosensor.

Fig.18 The 3D structure of ERα.

ERα contains six subunits, each of which has a LBD to accept the molecule of EE2. The LBD is folded into three-layered antiparallel α-helical sandwich comprising a central core layer of three helices sandwiched between two additional layers of helices. This helical arrangement creates a ‘wedge-shaped’ molecular scaffold that maintains a sizeable ligand-binding cavity at the narrower end of the domain. The remaining secondary structural elements, a small two-stranded antiparallel β-sheet and a helix are located at this ligand-binding portion of the molecule, and flank the main three-layered motif.

Fig.19 The 3D structure of ERα LBD.

Computational design of mutant ER receptors
In order to enhance the extent of the combination of ERα and EE2, position 421 was chose for virtual mutation. There are two important reasons to choose position 421. First of all, position 421 is an amino acid that exists as natural variants in different species, which should increase the probability that the resulting mutated protein is functional and properly folded. Secondly, the position 421 is relatively close to the ligand (just 5.9Å), and its mutation may alter the interaction and affect the binding ability.

Fig.20 The distance between ligand and the position 421.

Saturation mutagenesis was used with computational design. In order to evaluate of the effect of mutations, we used Ledock to calculate the binding energies, which brings much convenience to our model comparisons. LeDock is a simple proprietary molecular docking software that can be used for docking of ligands with protein target. Ledock adopts the algorithm of simulated anneal-genetic algorithm crossover to conduct conformational search. The docking scoring covers van der Waals interaction (VDW), electrostatic interaction (Coulombic interaction), hydrogen bond contribution and inter-molecular and intra-molecular clash of ligand. The detailed calculation is described as following.

Evaluation of binding free energy:

where,
∆E_ff--the interaction energy between the ligand and the protein calculated by the CHARMm force filed;
P_HB--hydrogen bonding penalty.

The equation for calculating ∆E_ff:

where,
∆E_vdW--the intermolecular van der Waals energy;
∆E_coul--the intermolecular Coulombic energy in vacuo;
∆E_strain--the strain energy of ligand upon binding;
∆G_solv--the change in solvation energy of ligand and protein upon binding.

The method for calculating P_HB:

where,
w--hydrogen bonding weights;

Table1 Hydrogen bonding weights

f_hb—the fraction of hydrogen bonding to that of an optimum geometry;

r--the distance between the hydrogen atom and the acceptor;
θ--the angle centered at hydrogen among donor, hydrogen and acceptor.
2.3 Result Discussion
The results about binding free energies with mutations are shown in Fig.4.

Fig.21 Evaluation of the position 421 mutations.

Simulation results show that mutants (M421C, M421T) may have better performances. In order to get a better understanding of mutations, Protein-Ligand Interaction Profiler(PLIP) was used which provides easy and fast identification of non-covalent interactions between biological macromolecules and their ligands. In addition, we want to screen mutants again with analysis of interactions.

Table2 Comparison of interactions between different mutants

Mutants	Hydrophobic Interactions		Hydrogen Bonds		π-Stacking		Interactions
	number	residue	number	residue	number	residue
Wide	5	350(ALA) 384(LEU) 387(LEU) 391(LEU) 525(LEU)	3	353(GLU) 394(ARG) 394(ARG)	1	404(PHE)
M421A	6	387(LEU) 391(LEU) 391(LEU) 404(PHE) 424(ILE) 424(ILE)	3	353(GLU) 394(ARG) 524(HIS)	1	404(PHE)
M421C	6	387(LEU) 391(LEU) 391(LEU) 404(PHB) 424(ILE) 424(ILE)	3	353(GLU) 394(ARG) 524(HIS)	1	404(PHE)
M421D	6	387(LEU) 391(LEU) 391(LEU) 424(ILE) 424(ILE) 525(LEU)	3	353(GLU) 394(ARG) 521(GLY)	1	404(PHE)
M421E	8	350(ALA) 384(LEU) 387(LEU) 391(LEU) 424(ILE) 424(ILE) 428(LEU)	3	353(GLU) 394(ARG) 521(GLY)	1	404(PHE)
M421F	9	350(ALA) 387(LEU) 391(LEU) 404(PHE) 421(PHE) 424(ILE) 425(PHE) 524(HIS) 525(LEU)	3	353(GLU) 394(ARG) 394(ARG)	1	404(PHE)
M421G	7	346(LEU) 387(LEU) 391(LEU) 391(LEU) 404(PHE) 424(ILE) 424(ILE)	3	353(GLU) 394(ARG) 524(HIS)	1	404(PHE)

As shown in the above mutants (Wide, M421A, M421C,M421D, M421E, M421F, M421G), analysis of all the mutants shows that most of the π-stacking remain unchanged. The hydrogen bond interaction was basically unchanged, but there were significant differences in the number of hydrophobic interactions among the mutants. In the interaction between ligand and protein, hydrophobic interaction plays a key role in stabilizing structure. Therefore, it is suggested to develop the M421F mutant which may enhance the degree of the combination between ERα and EE2.

Fig.22 PLIP(Wide).

Due to mutation, phenylalanine at position 421 can produce hydrophobic interaction with one end of the ligand, which can lift the one end of the ligand. Therefore, the ligand is closer to position 421 and closer to more amino acids, resulting in more hydrophobic interactions. More hydrophobic interactions make the ligand EE2 and the protein ERα binding stable.

Fig.23 PLIP(M421F).

In order to lower detection limit, we hope to achieve the desired detection performance by introduce by improving the combination between ligand and protein. Through the above analysis of LeDock and PLIP, using the strategy that mutating the position 421 M into F would show better detection results.

References

[1] Cowley SM, Parker MG. A comparison of transcriptional activation by ER alpha and ER beta. J Steroid Biochem Mol Biol. 1999 Apr-Jun;69(1-6):165-75. doi: 10.1016/s0960-0760(99)00055-2. PMID: 10418990.
[2] Brzozowski AM, Pike AC, Dauter Z, Hubbard RE, Bonn T, Engström O, Ohman L, Greene GL, Gustafsson JA, Carlquist M. Molecular basis of agonism and antagonism in the oestrogen receptor. Nature. 1997 Oct 16;389(6652):753-8. doi: 10.1038/39645. PMID: 9338790.
[3] Zhao H, Huang D (2011) Hydrogen Bonding Penalty upon Ligand Binding. PLoS ONE 6(6): e19923. https://doi.org/10.1371/journal.pone.0019923

To Top