Team:Melbourne/Model
Modelling
On this page we will be using a range of modelling techniques in order to attempt to gain a better understanding of the behaviour of the star peptide. Firstly, we shall aim to probabilistically determine the shape and structure of the star peptide. Using this stochastic model, we can estimate the diffusion time constant between arms of the star peptide via Brownian motion. Then, we will aim to compare the average diffusion time between star peptide arms with the average collision time of a normal solution and develop a relationship for the efficiency of the star peptide. This in turn will then be applied to the mevalonate pathway. Furthermore, as part of our collaboration with the National University of Singapore iGEM Team , we asked them to develop a simple diffusion simulator. This simulator will assist in gauging the behaviour of molecules surrounding the star peptides.
Characterisation of the Star Peptide
In this section we aim to develop a probabilistic model of the length of each arm of the star peptide, based on the number of amino acids in each arm. Obviously there will be complex interactions between the amino acids which will create a very unique protein structure. In order to simplify this protein folding problem, we will model the star peptide arms as self-avoiding random walks on a periodic lattice. There are a number of assumptions involved that will affect the accuracy of this model:
- Static Flexibility
- This refers to a random distribution of the bonds connecting amino acids. However, it is highly likely that neighbouring amino acids will sterically interact with one another to make the joining bond deterministic. Despite this, presence of these short-range correlations does not alter the basic random walk character of the chain statistics, it will simply affect the step size. This is because over a range of amino acids, the bonds will eventually become uncorrelated.
- Define
- $N_{seg}$: The number of amino acids in a segment, where bonds between segments are random
- $a$: The length of each segment
- $r_i$: The length of arm i, defined as the Euclidean distance from the centre of the star peptide to the end of the final amino acid of arm $i$
- Dynamic Flexibility
- The bonds between segments vary over time
- Self-avoiding random walk
- Self-avoiding means that once a point on the grid is occupied, it cannot be occupied again. It is effectively excluding a portion of the total volume.
MATLAB Simulation
A self-avoiding random walk simulation on a periodic lattice was constructed to visualise the structure of the star peptide. Using this simulation in MATLAB, we developed a probability distribution for r.
The star peptide shape model can be downloaded here.
The following values were used for the simulation:
- $N_{seg} = 100$
- $a = 1$
The simulation was run 10,000 times and the following probability distribution for the arm length $r$ was obtained.
The corresponding normal distribution curve is also plotted.
| Summary Statistics | |
|---|---|
| Mean | 12.57 |
| Standard Deviation | 2.43 |
The following is a visualisation of the star peptide's probabilistic diameter.
Theoretical Size
In order to obtain a theoretical relationship between arm length and the number of amino acids on each arm, we shall attempt to analyse the energy of the system. We can define the star peptide as four connected springs as shown. We will assume that each arm is the same and thus has the same spring constant.
Now, assuming the springs obey Hookean dynamics, we will develop a relationship between force applied and extension.
$$k_{eff} = (k+k)||(k+k)$$
$$k_{eff} = \frac{1}{\frac{1}{2k} + \frac{1}{2k}}$$
$$\therefore k_{eff} = k$$
From polymer physics, we know that:
$$F_{elastic} = kr$$ $$F_{el} = \frac{3k_B T}{Na^2} \cdot r$$ $$E_{el} = \int F_{el} dr$$ $$\therefore E_{el} = \frac{3k_B T}{2Na^2} r^2$$Repulsion from excluded volume (V),
$$E_{rep} = \frac{k_B T v N^2}{2r^3}$$Assume $v = a^3$,
$$E_{rep} = \frac{k_b T a^3 N^2}{2r^3}$$Calculating total energy,
$$E_{total} = E_{el} + E_{rep}$$ $$E_{total} = \frac{3k_B T}{2Na^2} r^2 + \frac{k_b T a^3 N^2}{2r^3}$$Finding the minimum $r$,
$$\frac{dE_{total}}{dt} = 0$$ $$\frac{3k_B T}{Na^2} r - 3\frac{k_b T a^3 N^2}{2r^4} = 0$$ $$\frac{r}{Na^2} = \frac{a^3 N^2}{2r^4}$$ $$r^5 = \frac{1}{2}a^5N^3$$ $$\therefore r = (\frac{1}{2})^{\frac{1}{5}}aN^{\frac{3}{5}}$$Substituting in values used in simulation ($N_{seg} = 100, a = 1$):
$$r = 2^{\frac{1}{5}}100^{\frac{3}{5}} \approx 13.8$$This theoretical value is slightly larger than the simulated value. For the rest of the model, we shall use the empirical value. Since we used $a = 1$, we can use the relationship $r = 12.57 \times a$.
Defining the Model
Our overall goal is to model the reaction kinetics of a solution where enzymes are connected to the star peptide and compare it to the kinetics of when the enzymes are floating freely in the solution. In the diagram below, the colour dots represent different enzymes. Assume:
- The concentrations of all enzymes and the star peptide are the same in both scenarios and are denoted by the constant $C$
- The substrates for the enzymes are floating around freely, and their concentration is a non-limiting factor
Scenario A
For a given enzyme concentration of $C$ molar, the mean square displacement (the average length between each enzyme) is:
$$ \begin{align} C M &= C \frac{mol}{L}\\ &= 1000C \frac{mol}{m^3}\\ &= 6.022 \times 10^{26} C \frac{particles}{m^3}\\ \end{align} $$ $$ \begin{align} MSD &= (\frac{1}{\frac{particles}{m^3}})^{\frac{1}{3}}\\ &= (\frac{1}{6.022\times 10^{26} C})^{\frac{1}{3}}\\ &= 1.18 \times 10^{-9} C^{-\frac{1}{3}} \end{align} $$Therefore, define the time constant $\tau_A$ as the average time of an enzyme substrate to diffuse between two enzymes in Scenario A:
$$\tau_A = \frac{MSD}{v} = \frac{1.18 \times 10^{-9} C^{-\frac{1}{3}}}{v}$$ Here, $v$ represents the velocity via Brownian motion given by: $$\frac{1}{2} mv^2 = \frac{3}{2}k_B T$$ $$\therefore v = \sqrt{\frac{3 k_B T}{m}}$$Scenario B
Assume that the steric interactions between the peptide arms interact so that the average star peptide shape is tetrahedral.
- $l$ = length between two enzymes in 3D
- $r$ = arm length
Therefore, assuming that all products from enzyme reactions continue flow towards the node where they are being consumed and don't interact with other star peptides, the length the substrates need to travel is $1.633r$. We can therefore define the time constant for Scenario B as:
$$\tau_B = \frac{1.633r}{v}$$Finding the reaction rate for Scenario B
As we go from Scenario A to Scenario B, the reaction rate will increase by a factor of $\frac{tau_A}{tau_B}$:
$$ \begin{align} \frac{tau_A}{tau_B} &= \frac{1.18 \times 10^{-9} C^{-\frac{1}{3}}}{v} \times \frac{v}{1.633r}\\ &= \frac{1.18 \times 10^{-9} C^{-\frac{1}{3}}}{1.633r}\\ &= 7.23 \times 10^{-10}(\frac{1}{rC^{\frac{1}{3}}}) \end{align} $$Mevalonate Pathway
We believe that a useful application for the star peptide is to increase the reaction kinetics of the mevalonate pathway.
This can be achieved with two star peptides joined together in conjunction with the enzymes attached to the ends of the arms as illustrated below.
Collaboration with the NUS Singapore iGEM Team
As part of our collaboration with NUS Singapore, we asked them to provide us with a diffusion simulator so that we could gain insight into the behaviour of molecules surrounding the star peptides.
The raw MATLAB files provided can be downloaded here.
The model was set up to mimic the two star peptides illustrated above. The aim of the model was to map the concentration fields of Acetoacetyl-CoA and HMG-CoA as they are produced and consumed at the arms of the star peptides. The following assumptions were used:
- the concentration of Acetyl-CoA is uniform and not a rate-limiting factor
- the initial concentration fields of Acetoacetyl-CoA and HMG-CoA are equal and uniform
- the diffusion constants of Acetoacetyl-CoA and HMG-CoA are equal
- the rate of conversion of Acetoacetyl-CoA to HMG-CoA at the HMG-CoA Synthase node is dependent upon the concentration of Acetoacetyl-CoA at the node
- the rate of consumption of HMG-CoA at the HMG-CoA Reductase node is dependent upon the concentration of HMG-CoA at the node
The adapted NU Singapore MATLAB files for this model can be downloaded here.
The following concentration fields were found:
These concentration fields clearly demonstrate the production and consumption of Acetoacetyl-CoA and HMG-CoA at the relevant nodes. However, the positioning of the contour lines brings into question a previously used assumption. The questionable assumption is that products from one enzyme reaction will flow towards the node of the other enzyme reaction where they are being consumed. If this were the case, then we would expect to see contour lines encapsulating and running between both the production and consumption nodes. Therefore, further testing using similar simulations should be performed in order to verify the accuracy of this assumption. Please note that the diamond shapes (which should be circles) are due to MATLAB's interpolation and not the diffusion simulator.
SimBiology Model
We shall now simulate the reaction kinetics of the mevalonate pathway using the SimBiology toolbox in MATLAB.
The SimBiology model can be downloaded here.
The mass action law was used to model the reaction kinetics of the normal cell. The mass action law multiplied by a factor of $\frac{\tau_A}{\tau_B} = 7.23 \times 10^{-10}(\frac{1}{rC^{\frac{1}{3}}})$ was used in the mevalonate pathway that uses the star peptide.
The following parameters were used:
- $C = 0.001M$
- $a = 4\unicode{x212B}$[1]
- $r = 12.57a = 50.28\unicode{x212B}$
The resulting reaction kinetics were:
