Introduction

Dynamic instability of microtubule during its lifetime is very interesting and worth researching. When a cancer patient is cured by taxol treatment fortunately, it should be owed to taxol's interaction with microtubule in dynamic growth process. So it's significant to explore principles of this process. We established a series of math works to clarify mechanisms of microtubule dilution system. Then we got data from wet parts of our project. Our model predicts observed data and phenomenon. In conclusion, we have presented a model that utilizes basic kinetics and it has been proved by experimental observations. We can use it to predict results in other conditions for guidance.

What we have done to analyze the dynamic process can be divided into three steps.

First we proved that taxol has a significant influence on microtubule growth by using the data which are provided from wet parts. Then we presented two differential equations model. One is an original differential equation model that can predict taxol’s influence on amount of microtubule. The other is a partial differential equation model expounding analytic solution. Since it’s hard to understand if someone lacks of necessary math skills, we build a program in MATLAB^@ to visualize the kinetics in microtubule disaggregation.

Modeling theory

Microtubule is made up of 13 protofilaments. Now there is a widely accepted feature about the microtubule that microtubule has highly complicated dynamic instability. Under vitro cultures conditions, on the one hand, subunits will polymerize automatically to form fine structure when subunits condition is above the critical concentration; on the other hand, the microtubule will depolymerize into subunits when subunits condition is under the critical concentration. Apart from that, the single microtubule will always in the stage of aggregation and disaggregation.

Tubulin is made up of two tubulin monomers which are nearly the same as each other. These two tubulin monomers are named α tubulin monomer and β tubulin monomer. Microtubule is made up of 13 protofilaments aggregated by tubulin dimers end to end. And microtubule can be the hollow tube with 13 protofilaments coiled into helix with each other, water in hollow part. The tube wall is 4～5nm thick. Tubulin dimers are incorporated into the growing lattice in the GTP-bound form and stochastically hydrolyze to GDP-tubulin, thus forming a GTP-cap. It is thought that the switching from growth to shrinkage occurs due to the loss of the GTP-cap.

Caplow M^[1] research shows that when the cap structure of microtubule plus end subunit containing GDP- beta tubulin instead of GTP- beta tubulin, microtubule becomes unstable and will quickly disaggregate.

As is shown in the figure, there are two kinds of Dimer, called GDP and GTP, with blue and red two connected to the circular. These dimers have close relationship with each other, and there are three important modes of their action.

1. GTP-tubulin dimer in endpoint can aggregate new GTP to make the single protofilament grow, and microtubules extend.
2. At the same time, the endpoint GTP may also be made off, thereby protofilaments shorter.
3. Any place of GTP (in addition to the right endpoints of the GTP) made made random hydrolyzed to GDP have a chance.

In order to verify and test the degree of the taxol’s influence on microtubule’s length.

We treated tubulin with taxol of different concentrations, measured their intensity of emission at around 520-530 nm wavelength.

First, we use the boxplot to conduct the descriptive statistical analysis. Then, we use the one-way analysis of variance to deal with the data given by the experiment.

1.0 - Descriptive statistical analysis

Before we use the one-way analysis of variance to deal with the data, we need to conduct some simple descriptive statistical analysis. This analysis can help us gain the overall understanding of the data, find the abnormal data, and then guide our next work.

We summarized tubulin-inflorescence fusion protein’s intensity of emission at around 520-530 nm wavelength under different taxol’s concentration treatments. The shape of experimental group emission curve and blank group emission curve are obviously different. That of experimental groups with different taxol’s concentration (final concentration: 0 μm, 0.5 μm, 5 μm, 25 μm, 50 μm, 100 μm, 250 μm and 500 μm) are basically similar except of the whole translation of wave band. Line graphs are shown here.

According to the picture above we can tell that adding taxol would influence the emission intensity of the fluorescent protein, which induces the shape and position changes of the intensity curve. Then we verified the influence of taxol to the fluorescent protein.

Because the interpreted variable is discrete, we use the Box Plot to describe the intensity of emission of microtubulin under the different concentration of the Taxol. According to the distribution of characteristic bands, we divide CPS value into four bands and use the Box Plot to describe them respectively. The specific four bands are as follows, 520nm-522.5nm, 522.5nm-525nm, 525nm-527.5nm and 527.5-530. And we get the following plot.

From the figure, we can see four boxplots. Each of them stands for a certain concentration of the TAX；the straight line in the box is the median of the data under a certain concentration, standing for the average of the light absorption value OD350 under this concentration. The upper line stands for the upper quartile of the data and the nether line stands for the nether quartile of the data. We can find that data have a big inner undulation under the concentration of TAX greater than 12.5, but when the concentration is between 5 and 12.5, the undulation of the data is tiny. In addition to this, the average of the light absorption value OD350 is high under the TAX’s concentration from 0 to 1 and the TAX’s concentration above 12.5.

2.0 - Theory of the one-way analysis of variance

By constructing the F-test statistics, we can use the one-way analysis of variance to study whether classification of the independent variable’s different levels can make significant influence on the variation of the continuous variable. If the levels have a significant influence, we can further give the 95% confidence interval of the dependent variable means under the different levels of the independent variable, and then we can analyze the degree of the different levels. But the precondition is that the data should satisfy the homogeneity of variance, in other words, the variance of the data should be the independent identically distributed. In the next part of the modeling, we will use the one-way analysis of variance to analyze the data, and then deal with the data.

To determine whether the data suitable for single factor analysis of variance

We use the SPSS to do the homogeneity test of variance with the data we got, the outcome is shown in the figure below:

From the figure, we can see the data’s variance is 0.519, bigger than 0.05. Therefore, we can think the data meets the requirement about the homogeneity of variance and we can use the one-way analysis of variance to deal with the data.

3.0 - Construct the F-test statistics

The independent variable is a classified variable which values 0 and 1 to describe whether the tax is added into the test tube. The dependent variable is the change of the micrutubule’ length, our modeling is shown below:

$$ y = u_i + \varepsilon_{ij} \tag{1} $$

In the equation (1), \(y\) is the dependent variable, the change of the microtubule's length. \(y_{ij}\) is the \(j\) observed value of the independent variable under the \(i\) level. is the mean of dependent variable under the \(i\) level. stands for the residual between dependent variable’s value and it’s mean value, also obey the normal distribution \(N(0, \sigma_i ^2)\)

Then we construct the F test statistics. First, we define the quadratic sum of the residual:

$$ SSE = \sum_{i=1}^k \sum_{j=1}^{n_i} (y_{ij}-\overline y_1)^2 \tag{2}$$

And the quadratic sum of the elements:

$$ SSA = \sum_{i=1}^k n_i (\overline y_{1}-\overline y)^2 \tag{3}$$

SSA reflects the variance between different levels and the difference is made by the different elements; SSE reflects the variance in a certain level and this random difference is due to the selected sample’s random. For example, the measured length of the microtubule will be different when we add taxol into the test tube.

On the basis of the theory, our F test statistics is:

$$ F = \frac{SSA/(n-k)}{SSE/(k-1)} \sim F(n-k, k-1) \tag{4}$$

The numerator of the equation is a part of the dependent variable which can be explained by the change of the independent variable, while the denominator of the equation can be explained by other random elements except the change of the independent variable. The proportion of the change of independent variable in all change of the dependent variable becomes bigger, in other words, F has a higher value, independent variable influence dependent variable more.

4.0 - The F-test on the data

Our F-test’s null hypothesis is:

$$ H_0:μ_1=μ_2=⋯=μ_k \quad vs \quad H_1=\quad not \quad H_0 \tag{5}$$

\(H_0\) stands for that different values of the independent variable（the added taxol‘s concentration) make no difference to the mean value of the dependent variable(the light absorption value OD350), in other words, the independent is not important to the dependent variable. Then we use R software to conduct F-test, the outcome is shown below:

From the outcome, \(F’s\) value is 1.337, nearly 1. Apart from that, the significance is 0.27, bigger than 0.05. Therefore, we can't think the independent variable make distinct influence on the dependent variable.

Differential Equation Model

For studying the dynamic progress of tubulin assembling under the influence of Taxol, we build the differential equations to describe such progress numerically.

First we have to consider what happened in the tubulin solution. Protein filaments will assemble into the tubulins automatically by arranging dimers. Specifically, the GTP-cap can aggregate or disaggregate at one side of the protein filament. Besides, there are chances that the GTP can hydrolysis to GDP, which can no longer combine the new GTP. In this process, the length and amount of the tubulin are in the dynamic equilibrium or changing.

We know that the differential equations can describe the numerical relationships and patterns between related promoters. Based on these studies, we can make the better plan, central or prediction in the experiments.

We consider the tubulin length is the function of time, assembling rate, disaggregation and hydrolysis rate. According to the literature, there are some relationships between these promoters.

Details of Original Differential Equation Model

a. - Model Assumption

1. We only consider tree main phenomenon in the solution: disaggregation, aggregation and hydrolysis.
2. There are three types of dimers, Taxol-dimer, GTP-dimer and GDP-dimer. Only Taxol-dimer and GTP-dimer, no matter they are dissociative or already fixed on the filament, can aggregate new Taxol-dimer and GTP-dimer.
3. Especially, the difference between Taxol-dimer and GTP-dimer is that new fixed aggregation Taxol-dimer will not disaggregate and become dissociative again.
4. According to literature ^[2], every Taxol-dimer and GTP-dimer may hydrolyze into GDP-dimer, which can no longer combine new dimers and growth.

b. - Parameter specification

c. - The explanation of the process

As we discussed in the modeling theory, we only consider three main phenomena: disaggregation, aggregation and hydrolysis.

Now we consider how these filaments growth specifically. At first, initial \(T_1\) and \(T_2\) are added in the solution. Then they aggregate into filaments forms, including \(G_{11}\), \(G_{21}\), \(G_{12}\), \(G_{22}\)(the meaning of these parameters are specified in table above). \(G_{11}\) and \(G_{21}\) can become \(G_{12}\) and \(G_{22}\) if new dimers are aggregated at the endpoint. Also \(G_{11}\)and \(G_{21}\) can disaggregate into the solution and become dissociative dimers (\(T_1\) and \(T_2\) again. If all the GTP-dimers at the endpoint become the GDP-dimers or there is not enough \(T_1\) and \(T_2\) in the solution, the growing process is stop.

The flow chart of this process can be illustrated as follow:

According to the flow chart, we can get the ODE equations as follow.

$\dot{T_1}=k_1 G_{11}-r_1(G_{11}+G_{21})$

$\dot{T_2}=-r_2(G_{21}+G_{11})$

$\dot{G_{11}}=-r_2(G_{11}+G_{21})+r_1(G_{11}+G_{21})-k_1 G_{11} \frac{G_{21}}{(G_{11}+G_{21})-h G_{11}}$

$\dot{G_{21}}=r_2(G_{11}+G_{21})-r_1(G_{11}+G_{21})+k_1 G_{11} \frac{G_{21}}{(G_{11}+G_{21})}$

$\dot{G_{22}}=r_2 G_{21}$

$\dot{G_{12}}=(r_1 + r_2)G_{11}-k_{1}G_{11}+r_1 G_{21}- h G_{12}$

$\dot{D}=h(G_{11}+G_{12})$

We can get numerical solutions easily by Matlab as soon as initial \(T_1\)and \(T_2\) are given. Obviously, initial \(T_1\) is related to the initial Taxol concentration. Accordingly, initial \(T_2\) is related to the initial GTP-dimer quantity. All these promoters are initial conditions in the experiment.

We consider a kind of initial conditions of \(T_1\) and \(T_2\). Actually, any initial \(T_1\) and \(T_2\) in the experiment can be modeled and we can predict the numerical trends. The solutions of one of our example are given in the following charts. (initial \(T_1\)=50000, initial \(T_2\)=3000, initial \(G_{11}\)=3000, Initial \(G_{21}\)=1000)

What we can get from our ODE model is \(G_1\) and \(G_2\), which indicate the quantity of filaments with GTP-cap and Taxol-cap, which is equal to the number of GTP-dimers at endpoint and the number of Taxol-dimers at endpoint. Sum of them is equal to the whole filaments units number in this dilution system. Since it should be proportionate to quantities of microtubules only by a coefficient \(\lambda \), we can use this to simulate the results whose condition is used in the wet parts.

Details of Partial Differential Equation Model

We know that the movement of any material is ruled by the certain laws of nature (physical and chemical laws). We establish some mathematical model including partial differential equations to describe these rules especially the mechanism of numerical relationships in the tubulin.

We build the partial differential equations and definition conditions based on the conservation laws.

The tubulin solution area \(\Omega\) is thought to be a kind of fluid motion. In the solution there are many assembling of tubulins, the hydrolysis, aggregations and disaggregation of the protein filaments. Although the mechanisms of these biological phenomena are still uncertain but we can build the equations based on conservation laws, which certainly rule the biomass.

We take a solution area as our study object. The quantity in the solution area is certainly ruled by conservation laws no matter what shape of the area is. Our models are based on conservation of mass.

We choose a solution area \(D\) in the \(\Omega\) . Considering in the time range , according to the mass conservation:

\(\overrightarrow{f_1}\) is the velocity of assembling of the solution which flow into the \(D\) . \(\overrightarrow{f_1}\) can be regarded as the product of the velocity \(\overrightarrow{v}\) and the assembling ratio \(k\) . However, there is no need to consider the \(\overrightarrow{v}\) and \(k\) . The inflow mass during \(dt\) through any \(dS\) of the \(\partial D\) can be given:

$$-\rho \cdot \overrightarrow{f_1} \cdot \overrightarrow{n} \cdot dS \cdot dt$$

\(\overrightarrow{n}\) is the unit normal vector toward outside so there is a minus in formula. \(\rho\) is the density of the tubulin solution in \(D\) . In a similar way \(\overrightarrow{f_2}\) is the velocity of the disaggregation mass which flow out the \(D\) and \(\overrightarrow{f_2}\) can be considered as the product of outflow velocity \(\overrightarrow{v_2}\) and disaggregation ratio \(r\) there is no need to consider the \(\overrightarrow{v_2}\) and \(r\) . So the outflow mass during \(dt\) through any \(dS\) of the \(\partial D\) can be given:

$$\rho \cdot \overrightarrow{f_2} \cdot \overrightarrow{n} \cdot dS \cdot dt $$

According to the mass conservation we can get the equation:

$$\iiint_D(\rho |_{t=t_2}- \rho |_{t=t_1})dxdydz=- \int_{t_1}^{t_2} dt \iint_{\partial D} \rho \overrightarrow{f_1} \overrightarrow{n}dS+ \int_{t_1}^{t_2} dt \iint_{\partial D} \rho \overrightarrow{f_2} \overrightarrow{n}dS$$

Assuming that \(\rho\) , \(\overrightarrow{f_1}\) , \(\overrightarrow{f_2}\) are all continuously differentiable, according to Gauss Formula we can get that:

$$\int_{t_1}^{t_2} dt [\iiint_D \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \overrightarrow{f_1})- \nabla \cdot (\rho \overrightarrow{f_2})]dxdydz=0 $$

Because of the continuity of integrand in and the randomicity of and , we can get that:

$$\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \overrightarrow{f_1})- \nabla \cdot (\rho \overrightarrow{f_2})=0, \quad \Omega \times (0, \infty) $$

And the definition conditions:

$$ \rho |_{t=0}=\rho_0$$

$$ \rho |_{t=t}=\rho_t $$

\(\rho_0\) , \(\rho_t\) can be measured directly in the experiment. For the \(\overrightarrow{f_1}\), \(\overrightarrow{f_2}\) , we can estimate the value range based on our experiments and determine the the most appropriate value according to the simulated annealing algorithm. Then the \(\rho(x,y,z,t)\) can be calculate according to the equation and definition conditions. Finally, the variation tendency can be described through the integration \iiint_D \rho(x,y,z,t)dxdydz

Visual Simulation

We applied to programing visualization in this complex process based on certain laws of Microtubule dynamic instability.

r=0.7 k=0.3 h=0.1

r=0.8 k=0.4 h=0.1