# Deriving kinetic equations for a complex reaction system

I have a system of protein reactions that I've derived differential equations for the kinetics of. I'm using them to run simulations to fit the rate parameters, but the simulations at present are working badly. I'd like to start by checking the system of differential equations. I realise they may be to complex for anyone to write them out here, but you know of software that could be used to check my work, please post it here. (A lot of people face similar problems)

The three monomers are $L$, $A_1$, $A_2$; the chemical equations for the monomers forming dimers are: $L + A_1 \rightleftarrows^{k_{a1}}_{k_{d1}} LA_1$, $\qquad$ $L + A_2 \rightleftarrows^{k_{a2}}_{k_{d2}} LA_2$,

and the dimers forming trimers are:

$LA_1 + L\rightleftarrows^{k_{a1}}_{k_{d1}} LA_1L$, $\qquad$ $LA_1 + A_2\rightleftarrows^{k_{a2}}_{k_{d2}} A_1LA_2$,

$LA_2 + L\rightleftarrows^{k_{a2}}_{k_{d2}} LA_2L$, $\qquad$ $LA_2 + A_1\rightleftarrows^{k_{a1}}_{k_{d1}} A_1LA_2$,

trimers forming tetramers:

$LA_1L+A_2\rightleftarrows^{k_{a2}}_{k_{d2}} LA_1LA_2$, $\qquad$ $LA_2L+A_1\rightleftarrows^{k_{a1}}_{k_{d1}} LA_2LA_1$,

$A_1LA_2 + L\rightleftarrows^{k_{a1}}_{k_{d1}} LA_1LA_2$, $\qquad$ $A_1LA_2 + L\rightleftarrows^{k_{a2}}_{k_{d2}} LA_2LA_1$,

dimers forming tetramers

$LA_1+LA_2\rightleftarrows^{k_{a1}}_{k_{d1}} LA_1LA_2$, $\qquad$ $LA_1+LA_2\rightleftarrows^{k_{a2}}_{k_{d2}} LA_2LA_1$,

tetramer molecules changing their configuration

$LA_1LA_2 \rightleftarrows^{k_{a4}}_{k_{d4}} B$, $\qquad$ $LA_1LA_2 \rightleftarrows^{k_{a3}}_{k_{d3}} B$.

Could anyone either point me in the direction of software to derive the system of differential equations, or do the derivation?

• I'm afraid you should know how to do this yourself. Also if you want to only integrate ODE's you shouldn't worry about CPU time. May 2 '16 at 22:30
• I'd use a matrix approach. Least likely to make errors. Set up the stoichiometric matrix. May 3 '16 at 9:55
• I'm assuming that you want time profiles not just equilibrium values. In this case use a Master Equation approach to set up a matrix and solve for eigenvalue values & eigenvectors numerically. Use these to reconstruct the time profiles etc. You can use Maple or Mathcad, both commercial (& expensive), or use Python which is free, but just as easy to use for numerical work. Oct 22 '16 at 11:58

In many cases, complex systems of even simple differential equations do not have simple and obvious solutions. Therefore, it could be too difficult and not reasonable to find fully analytical solutions.

Usually, numerical solutions can be easily calculated using iterative algorithms. However, such computer codes can be susceptible to calculation artefacts. Furthermore, the validation and verification of the model and the code can be a difficult issue.

I made a compromise and used an analytical approach when solving systems of differential equations, but I allowed the code to evaluate subexpressions involving constants numerically. Thus, I made numerical calculations to combine all constant parameter values as soon as possible in order to limit the number of parameters in the equations (using excessive digits in order to avoid accumulation of round-off errors).

I found that this approach can be easily and quickly implemented using the commercial software Maple. However, other similar mathematical software could be suitable, too.

I successfully used this approach for the evaluation of radioactive decay of the inventory of a nuclear reactor core in dependence of time after shutdown. The calculation of the radioactive decay of the relevant radionuclides involved more than 100 individual differential equations, which mainly describe chains of simple first-order reactions, e.g.:

$$\frac{\mathrm da_{\text{I-131}}(t)}{\mathrm dt}=-\lambda_{\text{I-131}}\cdot a_{\text{I-131}}(t)+0.778\cdot\lambda_{\text{I-131}}\cdot a_{\text{Te-131m}}(t)+1\cdot\lambda_{\text{I-131}}\cdot a_{\text{Te-131}}(t)$$

> de_a_I131:=diff(a_I131(t),t)=-l_I131*a_I131(t)+0.778*l_I131*a_Te131m(t)+1*l_I131*a_Te131(t);


A large number of further differential equations could be added in order to model transport phenomena during postulated accidents (e.g. loss of coolant).

Differential equations that were dependent on each other and the corresponding initial conditions were grouped into systems of differential equations, e.g.:

> sys_131:={de_a_Xe131m,de_a_I131,de_a_Te131,de_a_Te131m,a_Xe131m(t0)=core_Xe131m,a_I131(t0)=core_I131,a_Te131(t0)=core_Te131,a_Te131m(t0)=core_Te131m};


The systems of differential equations were solved analytically and the subexpressions involving constant parameter values were evaluated numerically.

> sol_131:=evalf(value(dsolve(sys_131,[a_Xe131m(t),a_I131(t),a_Te131(t),a_Te131m(t)])));


The results could be easily rearranged into a suitable form, e.g.

\begin{aligned}a_\text{I-131}(t) ={} & -1.6030637832147524320\cdot10^{17}\cdot\mathrm e^{-6.4180294496291232354\cdot10^{-6}\cdot t}\\&-7.9769363886868464215\cdot10^{15}\cdot\mathrm e^{-4.6209812037329687295\cdot10^{-4}\cdot t}\\&+4.7213833147101620896\cdot10^{18}\cdot\mathrm e^{-9.9782796169607015473\cdot10^{-7}\cdot t}\end{aligned}

The results could also be directly used for further calculations within the same worksheet.