Team:SYSU-Medicine/Injection System
<!doctype html>
Modeling Part 2
Injection System
In this section, our experiment focused on the medicine effect on the mouse. How can we measure the complex relationship between actinomycin D and M1? How does the concentration of nitroreductase grow with time? What is the best time and best prodrug volume for injection? Our models tried to solve these problems.
The process of our treatment is complex. Because it contains the microcosmic process, in which M1, uninfected and infected tumor cells interact with each other. Besides it is related with macroscopic process, in which the prodrug quantity of injection has a complex relationship with the initial concentration of prodrug in partial area. So we separate the whole process into partial process and injection process, and then we analyze them detailedly.
Partial Process
Tumor cells often gather in the specific area, that is why we choose the partial process for our research.
In the partial process, the initial condition is the specific concentration of M1 and actinomycin D
prodrug as the pro-medicine, along with the initial concentration of tumer cells. So we have the
initial condition
U(0) = φU, I(0) = 0, M(0) = φM, cm(0) = φm, cE(0) = 0
Then M1 will infect the tumor cells specifically. While in the tumor cells, M1 not only makes self replication but also translate the nitroreductase. The nitroreductase is so big that it can not be transported outsides. Luckily, our prodrug can enter the tumor cells due to its special perssad. So prodrug can be restored to actinomycin D under the help of nitroreductase. Actinomycin D producted could damage the tumor cells and at the same time, boost the production of nitroreductase.
Tumor And M1
The process of tumor cells with M1 is similar to the process of veros with M1, so we obtain the ODE
Model in the section 1. But we improve the model considering the effect of actinomycin D:
$$\left\{ \begin{array}{lcl}
\frac{dU(t)}{dt}=rU(t)(1-\frac{U(t)+I(t)}{K})-\beta U(t)M(t)-k_1 c_m(t)U(t)\\
\frac{dI(t)}{dt}=\beta U(t)M(t)-\alpha I(t)-k_2c_m(t)I(t)\\ \frac{dM(t)}{dt}=\sigma
I(t)-uM(t)+\delta c_m(t) \end{array} \right.$$
In the improved ODE model, k1 is the coefficient of actinomycin D killing uninfected veros, while k2 is the coefficient of actinomycin D killing infected veros. In addition, let δ be the coefficient of actinomycin D boosting the growth of M1. The other coefficients are the same as in Section 1.
Traditional Michaelis-Menten Equation
Now for the harder part. Our aim is to find out the influence of the changing concentration of the nitroreductase to the partial process. To find out the relationship of enzyme and reactant, we obtain the Michaelis-Menten Equation. In our experiment, the enzymatic reaction condition is suitable for the Michaelis-Menten Equation, especially for the assumption of Briggs-Haldane improved model.
E + S ↔︎ ES → E + P
So there is convincing to use the Michaelis-Menten Model.
Let cp be the concentration of prodrug,
while cm is the concentration of
actinomycin D. cp and cm are the total concentration ctotal,
ctotal = cm + cp
for the ratio of reactant and product in our chemical equation is one to one. First we assume that
the concentration of nitroreductase is constant, so the model is
$$\frac{dc_p}{dt}=-\frac{dc_m}{dt}=-\frac{\beta_1
c_p}{\beta_2+c_p}$$
.
The integral of the equation is below (we have considered the initial condition of our model).
$$\left\{ \begin{array}{lcl} c^{\beta_2}_p e^{c_p}=e^{-\beta_1(t+c^*)}\\
(c_{total}-c_m)^{\beta_2}e^{c_{total}-c_m}=e^{-\beta_1 (t+c^*)}\\ c^*=-\frac{\beta_2 ln
c_{total}+c_{total}}{\beta_1} \end{array} \right .$$
then we get the figure of
traditional Michaelis-Menten Equation that
Improved Michaelis-Menten Equation
However, the concentration of nitroreductase is time-varying because it can be produced by the
M1-infected tumor cells. Based on the usual method of Michaelis-Menten Equation and for the
convenience as well, we assume that the coefficients of Traditional Michaelis-Menten Equation β1 and β2 are the affine transformation of the
concentration of nitroreductase cE. So
the differential form of our improved Michaelis-Menten Equation is
$$\frac{dc_p}{dt}=-\frac{dc_m}{dt}=-\frac{(\beta_1+\gamma_1 c_E)
c_p}{(\beta_2+\gamma_2 c_E)+c_p} .$$
Then we obtain the integral form of the improved Michaelis-Menten Equation as
$$\left\{ \begin{array}{lcl} c^{\beta_2+\gamma_2 c_E}_p
e^{c_p}=e^{-(\beta_1+\gamma_1 c_E)(t+c^*)}\\ (c_{total}-c_m)^{\beta_2+\gamma_2
c_E}e^{c_{total}-c_m}=e^{(-\beta_1+\gamma_1 c_E) (t+c^*)}\\ c^*=-\frac{(\beta_2+\gamma_2 c_E) ln
c_{total}+c_{total}}{\beta_1+\gamma_1 c_E} \end{array} \right .$$
The figure of the equation is below.
Integrated Partial Process ODE Model
Now, we have to consider the influence of the time-varying concentration of nitroreductase to the process.
We know that the concentration of nitroreductase can rise up just due to the M1 infection. And it will decrease because of the death of infected tumor cells. These are two main ways to influence the concentration of nitroreductase.
Finally, we get the ODE model of partial process as follows.
$$\left\{ \begin{array}{lcl}
\frac{dU(t)}{dt}=rU(t)(1-\frac{U(t)+I(t)}{K})-\beta U(t)M(t)-k_1 c_m(t)U(t)\\
\frac{dI(t)}{dt}=\beta U(t)M(t)-\alpha I(t)-k_2c_m(t)I(t)\\ \frac{dM(t)}{dt}=\sigma
I(t)-uM(t)+\delta c_m(t)\\ \frac{dc_m}{dt}=\frac{(\beta_1+\gamma_1 c_E)
(c_{total}-c_m)}{(\beta_2+\gamma_2 c_E)+ (c_{total}-c_m)}\\
\frac{dc_E(t)}{dt}=k_3M(t)-k_4(\alpha I(t)+k_2 c_m(t)I(t)) \end{array} \right .$$
In our equation group, the first equation ends with the element containing U(t) to avoid the situation that the rate is under
zero. It is the same in the second equation. The initial condition of our equation group is
U(0) = φU, I(0) = 0, M(0) = φM, cm(0) = φm, cE(0) = 0.
Thus, we obtain the ode45 or ode15s to solve the equation group numerically.
Parameter Estimation
The number of our equation group is big in some way. Some of the parameters can be estimated by other methods (such as measurement), but some of them can only be estimated by fitting the data.
In the Improved Michaelis-Menten Equation, β1 and β2 are related with the concentration of nitroreductase linearly. So we use the linear regression model to fit the data. We should first measure the 5 min rate under different concentration of nitroreductase and then fit the data with traditional model in section 2.1.2, in order to get the estimated parameters β1 and β2 that exist in the traditional Michaelis-Menten Equation. Next we use linear regression model to find out the relationship between cE and β1, β2. Besides our method, there are Lineweaver-Buck Method and Wilkinson Method. However, the error in our method is the lowest among three.
In our experiment, the prodrug production is low. We haven’t measured the concentration of drug under different concentration of nitroreductase. However, we will finish the parameters estimation of our model based on the method we used above. We analyse the model based on different paramter values. We first analyse the effect of β1, and find that it only affects the maximum of the curve. β2 affects the rate of getting to the maximum.
Injection Process
In this section, we will take the process of prodrug from injection to targeting area into consideration. In addition, we only consider the situation in immune-free mouse in the one-injection curve.
Indicator
To measure the degree of pesticide effect, kill rate is an indicator that is often used in the medicine research. In our research, instantaneous kill rate is more important, for it can change greatly, catching our attention. So we all consider instantaneous kill rate as our index. We assume that kill rate is related with the concentration of M1 and actinomycin D. Theoretically, M1 and actinomycin D kill tumor cells independently. However, they influence each other actually. So we have to add the method of Elastic Net or interation term to eliminate the autocorrelation of them (we choose the latter one). Furthermore, kill rate has positive correlation with the concentration of tumor cells (in our model we use ratio of survived tumor cells), so we finally show that
R̂kill = (β1cactin + β2cM1 + β3cactincM1 + β4)Rtumorα, 0 < α < 1
To guarantee value of kill rate in [0, 1], we use the function ϕ(x) = 1 − e − x. In our experiment, the parameters are estimated by the data in confirming the validity of actinomycin D boosting the therapy of M1. The parameters are β1 = 0.6615, β2 = 0.3418, β3 = 1.9844, β4 = 0 and α = 1. Then we get the kill rate function and the figure.
Rkill = 1 − e − (0.6615cactinomycin D + 0.3418cM1 + 1.9844cactinomycinDcM1)Rcancer
Compartment Model
To estimate the time-varying distribution of prodrug, we obtain the compartment model. A very important assumption in our compartment model is that medicine follows the uniform distribution in a room.
Thus we also use the ODE model according to the relationship in the process.
$$\left\{ \begin{array}{lcl}
\dot{x}_1(t)=-k_{12}x_1-k_{13}x_1+k_{21}x_2+f_0(t)\\ \dot{x}_2(t)=k_{12}x_1-k_{21}x_2
\end{array} \right.$$
We know that the concentration c(t) multipling the
volume V is the quantity x(t).
xi(t) = Vici(t), i = 1, 2
And then we can get that
$$\left \{ \begin{array}{lcl}
\dot{c}_1(t)=-(k_{12}+k_{13})c_1+\frac{V_2}{V_1}k_{21}c_2+\frac{f_0(t)}{V_1}\\
\dot{c}_2(t)=\frac{V_1}{V_2}k_{12}c_1-k_{21}c_2 \end{array} \right.$$
In our experiment, the type of injection is rapid intravenous injection. So f0(t) and the initial condition is
$$\left \{ \begin{array}{lcl} f_0(t)=0\\ c_1(0)=\frac{D_0}{V_1}\\
c_2(0)=0 \end{array} \right .$$
Then we could get the solution of our equation that
$$\left \{ \begin{array}{lcl} c_1(t)=Ae^{-\alpha t}+Be^{-\beta t}\\
c_2(t)=\frac{D_0 k_12}{V_2 (\beta-\alpha)}(e^{-\alpha t}-e^{\beta}) \end{array} \right
.$$
The parameters in the model are
$$\left \{ \begin{array}{lcl}
A=\frac{D_0(k_{21}-\alpha)}{V_1(\beta-\alpha)}\\ B=\frac{D_0(\beta-k_{21})}{V_1(\beta-\alpha)}\\
\alpha+\beta=k_{12}+k_{21}+k_{13}\\ \alpha \beta=k_{21}k_{13} \end{array} \right
.$$
Then we can get the time-varying concentration of total medicine.
Parameter Estimation
In our experiment, the concentration of actinomycin D is measured under different time. We can see the data in the table.
| ti | 1 | 2 | 4 | 8 | 16 | 24 | 36 |
|---|---|---|---|---|---|---|---|
| c(ti)μg/mL | 0.1005 | 0.0845 | 0.0638 | 0.0349 | 0.0136 | 0.0044 | 0.0020 |
Then we use the former 4 iterations and latter 4 iterations to suit the model ln c(t) = ln P1 − P2t and we finally get that A = 0.09154, α = 0.1204, B= and β = − 0.009719. And we calculate other parameters that are shown in the coming table.
| A | 0.09154 | k21 | 0.1337 |
|---|---|---|---|
| α | 0.1204 | k13 | 0.1356 |
| B | 0.1159 | k12 | 0.0017 |
| β | 0.1506 |
The figure that contains the nonlinear fitting method is below.
Then we could see that the concentration of mixed prodrug changes as time goes by. If we only inject the mouse once, the time-varying concentration of mixed prodrug and Half-Life Period are provided in the figure. We also provide the concentration of M1 both in the serum and tumor area. Compared with the data, our model curve is suitable for our research.
Integrated Injection Model
Now we have the Partial Process and Injection Process, and then we want to mix them into one process. The concentration of total medicine in the tumor area is only related to the Compartment Model, and limits the time-varying partial process. Thus we could get the curve of concentration of tumor cells, titer of M1, concentration of nitroreductase, and concentration of actinomycin D in blood system and in tumor area as follows.
As for the instantaneous kill rate, we first analyse the data we got. Here is the relationship matrix of the concentration of uninfected tumor cells, infected tumor cells, titer of M1 and the concentration of actinomycin D. The data is from our model. Searching for the relationship between two variable just based on our model can not be convincing. However, there is rule between them which can be confirmed not just from references but from our experiment. Thus, we consider the possible error in our model by adding Gaussian noise.
We can see that there are strong relationship between variables. The relationship above confirm our assumption. Thus, the index Instantaneous Kill Rate is reasonable. And we get the time-varying Instantaneous Kill Rate as follow.
Finally, we could get the injection regime by adapting the curve to the upper and lower boundary. After trying a lot of times, we find that 0.58μg/mL is the best injection quantity in the experiment. (The value of best injection time is a little different in this figure compared to the above figure, because we use the Runge Kutta Algorithm to solve the impulsive differential equation, which can be less complex than the exact solution in coding.)
Video of 2-dimentional Distribution
To see the time-varying distribution of M1 inside mouse, we select 6 photos of 2 mice with fluorescent and a photo of these 2 mice without fluorescent. We convert every image into an RGB matrix. Then we try to transform the 6 matrixs and reduce the target dimension. Next we try to use BP neural network to fit the time-varying image and create 121 photos. Finally we make these images into movie.
Video of 2-dimentional Distribution
Renference
[1]Jiang Qiyuan, Xie Jinxing, Ye Jun (2010). Mathematical model (4th ed.). Beijing: China Higher Education Press, 338-346
[2]Yang bowen, liu ping. Pharmacokinetics two-compartment model [J]. Journal of natural science of Harbin normal university,2016,32(01):13-15+66.
[3]Yuan jin, Shi Lei, Zhao Shujin. Pharmacokinetics and Compartment Model Parameters of the Two-Compartment Model for External Vascular Administration Based on Excel Function [J]. China Pharmacy,2008(02):106-108.
[4]https://2016.igem.org/Team:SYSU-MEDICINE/Modeling