Team:DUT China A/Model
Introduction
1. What is our model
We use mathematical modeling to determine parameters and narrow our range of options. Specifically, three models were adopted for our project: The binding of aptamer to target sites，Determination of inhibitor and the appropriate time for growth of hydrogels.In the first model，we establish a differential equation based on the molecular mechanism in cell micro-environment and solving them.
In the second model, we interpreted the issue as an optimization problem and used thermodynamics theory to quantify objective function and constraint function.
In the third model, we use post feedback neural network to predict the radius of colloidal particles in given time. By training the network with experimental data, we can get an approximate mapping from time to radius of colloidal particles.
2. Why modeling?
The main purpose of our modeling is to determine or narrow some options for experiments. For example, the determination of inhibitor, the binding of aptamer to target sites, and the appropriate time for growth of hydrogels need to be modeled. With the preliminary theoretical reference, the members in the experimental group can do the experiment effectively.3. What help does the experiment provide for the model?
The experiments provide us with necessary data, such as MFI in the binding of aptamer to target sites and the data of micelle radius in the appropriate time for growth of hydrogels. The former helps us determine the parameters and verify the model, and the latter helps us train the network.4. How does the model guide the experiment？
The process of solving the model parameters guides the experiment. For example, we need to determine the reaction rate constant in the binding of aptamer to target sites. Therefore, we did the experiment of target binding to aptamer by fluorescence labeling on aptamer. The residual concentration of the aptamers in the reaction equilibrium is obtained by fluorescence intensity, and the rate constants can be solved by it.The binding of aptamer to target sites
In our experiments, an important step is to combine the aptamers with the cell targets, and subsequently produce hydrogels. However, Cells will not be fully reacted if they react for too short time and be damaged if they react for too long time, which both affect the formation of hydrogels. Therefore, it is very important to find a suitable combination time.
The binding reaction between aptamer and target spots is the first-order kinetics：
aptamer+target↔a-t
There is a chemical positive reaction rate equation: $$\begin{equation} v_1=kc_a(t)c_t(t) \end{equation} $$
(2.1)
Similarly, the chemical inverse reaction can be described by the following one.
$$\begin{equation}
v_2=\frac{1}{k}c_{a-t}(t)
\end{equation}
$$(2.2)
\(c_a(t)\), \(c_t(t)\) is the residual concentration of aptamer and target over time t. k is the constant of chemical reaction rate. Therefore, it can be considered that the total reaction rate v=\(v_1\)-\(v_2\). When v = 0, the chemical reaction is in equilibrium. Then we can get a differential equation:
$$\begin{equation}
-\frac{dc_a(t)}{dt}=-\frac{dc_t(t)}{dt}=v=v_1-v_2=kc_a(t)c_t(t)-\frac{1}{k}c_{a-t}(t)
\end{equation}
$$(2.3)
Through the relationship among \(c_t (t)\) and \(c_(a-t) (t)\), we know that:
$$\begin{equation}
c_a(t)-c_a(0)=c_t(t)-c_t(0)=-c_{a-t}(t)
\end{equation}
$$(2.4)
Then, based on the equation (1.1), (1.2), (1.3) and (1.4), we can get the following equation：
$$\begin{equation}
-\frac{dc_a(t)}{dt}=k[c_t(0)-c_a(t)+c_a(t)]c_a(t)-\frac{1}{k}[c_a(0)-c_a(t)]
\end{equation}
$$(2.5)
The solution is：
$$\begin{equation}
t=c_0-\frac{1}{k(x_1-x_2)}ln(\frac{c-x_1}{c-x_2})
\end{equation}
$$(2.6)
\(x_1\) and \(x_2\) are two roots of equation \(\begin{equation}
Kx^2+(mk+\frac{1}{k})x-\frac{c-a(0)}{k}=0
\end{equation}
\). \(\begin{equation}
m=c_t(0)-c_a(0)
\end{equation}
\), can be calculated from the initial concentration. Substituting the initial value into the equation (2.6)，the constant can be calculated: \(\begin{equation}
c_0=\frac{1}{k(x_1-x_2)}ln(\frac{c_a(0)-x_1}{c_a(0)-x_2})
\end{equation}
\)In the experiment, we added the concentration of target point, \(\begin{equation} c_t(0)=8.3*10^{-10}mol/L - 1.6*10^{-8}mol/L \end{equation}\), which had little effect on M, so we finally confirm \(\begin{equation} c_t(0)=10^{-8}mol/L \end{equation} \), \(\begin{equation} c_a(0)=2*10^{-7}mol/L \end{equation}\). According to the experimental data, when the reaction is in equilibrium, \(\begin{equation} c_a=115*10^{-9}mol/L \end{equation}\).
When the reaction is in equilibrium，we can consider the time as infinity.
(2.7)
The solution is：$$k=1.5567* 10 ^3$$
Based on Table 1 we calculated the fitted curve between concentration and MFI.y=16.05x+110.31.y is concentration, while x id MFI. According to the function and the experimental data of MFI over time, we compute the changes of concentration over time. The specific results are shown in the Table 2.
Table 1: concentration with MFI
Table 2: concentration over time
By substituting K into the equation (1.6), the following results can be obtained (And compared with the experimental results):
Figure 1: The change of aptamer concentration with time
It can be seen from the figure that the experimental data and the theoretical simulation value are basically consistent. According to the theoretical calculation value, when the reaction time reaches 1 hour, the reaction balance is basically reached. We can take one hour as the best time of experiment.
Conclusion
We set up an ordinary differential equation with parameters through the reaction rate equation, and solved the parameters through the experimental data, and determined a suitable reaction time of 1 hour for the experiment. This model can be applied to other first-order kinetics, the parameters can be determined by similar methods, then the theoretical value of the reaction can be obtained.Determination of inhibitor
The reaction of aptamer binding to target is also the reaction of primer exposure. The exposure of Primer is directly related to the formation of hydrogel, so we need a suitable cover length to make the Primer stable enough before aptamer binds to the target, and to make the Primer smoothly exposed when aptamer binds to the target. In this way, we can compute a feasible region where we can optimize.
3.1 Objective function
Aptamer is a single strand of DNA that is easy to have the phenomenon of self-rotation. In order to suppress the self-rotation of DNA, we need to maximize the length of the masking sequence in the feasible domain.2, The establishment of feasible domain
The assumption of nearest neighbor method [1] :
The initial concentration of two single nucleic acid chains in the hybridization reaction is the same;
∙∙The nucleic acid hybridization reaction satisfies the distate model, i.e., the dechain and the formation of the double strand are always dynamically balanced;
∙∙The hybrid nucleic acid double strand is a single double helix structure.
3.2 The first of the constraint function
Based on the assumption, we can compute the model in this way. This reaction can be divided into two different reaction. ①The first reaction is between aptamer and target protein: $$\begin{equation} AT \leftrightarrow A+T \end{equation} $$ The equilibrium constant of the reaction is:$$\begin{equation} K_d1=\frac{c_ac_t}{c_at} \end{equation} $$ ②Another reaction is between aptamer and inhibitor $$\begin{equation} AI \leftrightarrow A+I \end{equation} $$ and the equilibrium constant of the reaction is: $$\begin{equation} K_d2=\frac{c_ac_I}{c_at} \end{equation} $$ The whole reaction is: $$\begin{equation} AC+P \leftrightarrow AP+C \end{equation} $$ To describe the direction and stability of the reaction, we introduce the thermodynamic concept of free energy. $$\begin{equation} \ ΔG_c=-RTln(K) \end{equation} $$
(3.1)
\( ΔG_c\) is the free energy of the enzyme, R is the thermodynamic constant, T is the temperature, and K is the equilibrium constant of the chemical reaction. \(\begin{equation}
K=\frac{K_d1}{K_d2}
\end{equation}
\)
When the free energy of aptamer self folding is smaller than the complementary energy of aptamer and masking sequence, aptamer self folding can be prevented. So the first constraint function is:
$$\begin{equation}
\sum_{i=1}^{k-1}\ ΔG_i\leq\ ΔG_c
\end{equation}
$$(3.2)
\(ΔG_i\) is the \(ΔG\) of i-th and (i+1)-th base-pairs’ position. The different kinds of base-pairs’ Gibbs energy are shown below. K is the length of the inhibitor.
Figure2：Base-pairs’ Gibbs energy
3.3 The second constraint function
According to the thermodynamic model, in order to expose the inhibitor, it’s necessary that the Gibbs energy of aptamer is greater than that of inhibitor. The compute of the energy is another constraint condition.$$\begin{equation} \sum_{i=1}^{k-1} ΔG_i\geq\ ΔG_a \end{equation} $$
(3.3)
\(ΔG_a\)=-6.25, is the Gibbs energy of aptamer, computed by software oligo analyzer.
3.4 Solution
The results show that the inhibitor length of aptamer is 5-10. Considering the multi-objective optimization, and trying to prevent DNA self-folding, we found that the suitable number of inhibitor are 8 and 9 by experiment, so we finally chose 9 bases as the inhibitor of aptamer part, and the inhibitor of linker part is not limited, so we took 6 bases. Then we determined the masking sequence as TCAGGC-CACTACAGABecause the second constraint function needs to use the oligo analyzer,we just show a built-in code to calculate the first constraint. Just enter the sequence in uppercase.
3.5 Conclusion
The assumption of nearest neighbor method provides the theoretical basis for the quantification of the constraint function. It is the establishment of constraint functions that allow us to transform it into an optimization problem. Because of this model’s small computation, it’s easy to confirm the best length by exhaustive method.The basic idea is from another IGEM team in Xiamen university. Based on the existing results, we add a restriction function to prevent DNA single strand from self folding and find out the feasible region of inhibitor. Finally, we determined the inhibitor by experiments as follows: TCAGGC-CACTACAGA.
The appropriate time for growth of hydrogels
Hydrogel is a large number of DNA single chains generated by DNA rolling ring amplification. Because of its good structural characteristics, we use it as prison clothes. Because the whole environment is not suitable for cell survival, too long time will lead to cell death, and too short time will lead to failure to form, it is important to determine the forming time of hydrogel. First, we consider establishing differential equations to solve this problem, but the rate of DNA rolling ring amplification is difficult to express (There are many effects that are difficult to quantify, such as DNA degradation), and in order to describe this process theoretically, we still describe the derivation of differential equations. Finally, we decided to adopt the BP neural network to approximate the mapping from time to the radius of the hydrogel to solve it.
4.1 Establishment and solution of network
We normalize the data and use two hidden layers in network, the number of neurons in each layer is 3, the active function from input layer to hidden layer is logsig, and the active function from hidden layer to output layer is tansig. The specific structure is as follows:
After training with experimental data (Table 3), the mapping of the network is:
Figure 3: The mapping of the network
So we get the suitable time:8 hours - 8 hours and 20 minutes.
Table 3: The radius of hydrogels particles over time
4.2 Derivation of differential equation
Notion:A and B are linear structures of complementary ring chains
\(n_1\): The number of whole circle sequence in DNA chain A
\(n_2\): The number of whole circle sequence in DNA chain B
\(c_1\): The number of DNA chain A
\(c_2\): The number of DNA chain B
\(V(n_1,n_2,t)\): The rate of DNA rolling ring amplification
Considering the simultaneous reaction of multiple chains, the growth rate of the number of rings of A(B) is directly proportional to the number of base points, the ratio is V, where V is the reaction rate of a single chain, which is related to N1, N2 and time. The basic point is the starting position of the chain, and the number of base points of A and B is n_2-c_2+1, n_1-c_1.Similarly, the growth rate of A (B) chain number is directly proportional to B (A) chain number, and the ratio is V[2][3].
So we can get:
4.3 Conclusion
We use experimental data and BP neural network to fit the mapping from time to hydrogel particle diameter. The suitable time for rolling ring amplification is 8 hours-8 hours and 20 minutes. After that, we set up the differential equations. Theoretically, knowing the change of the rate V about N1, N2 and T, we can numerically solve the change of the number of bases with time, but we did not find this relationship. We hope that our later research can make a breakthrough in this part.Reference
[1] JOHN SANTALUCIA, JR*, A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics, Proc Natl Acad Sci Anal Chem, 1998, 95, 1460-1465.
[2] Gevertz, J. L., Dunn, S. M., & Roth, C. M. (2005). Mathematical model of real-time PCR kinetics. Biotechnology and Bioengineering, 92(3), 346–355. doi:10.1002/bit.20617
[3] Beals T P , Smith J H , Nietupski R M , et al. A mechanism for ramified rolling circle amplification[J]. BMC Molecular Biology, 2010, 11(1):94.