ODE

ODEs applicability to enzyme kinetics

Ordinary Differential Equations (ODEs) are at the heart of modeling enzyme kinetics. They are equations describing the rate of change of a variable e.g. product concentration. To illustrate their applicability consider the following reaction scheme:

Consider the rate of change of the enzyme complex (AE_A). There are 3 possibilities:

$$A+E_A \rightarrow AE_A $$ $$AE_A \rightarrow A + E_A$$ $$AE_A \rightarrow H + E_A$$

It can break down into H and E_A or it could break down into A and E_A or it could be formed by A and E_A. The first two possibilities will happen proportional to the amount of AE_A. the latter possibility is proportional to the product of the concentrations as it’s proportional to the number of collisions.

Assigning rates to these steps which represent the probability of each process happening you would therefore get the following equation: $$\dot{AE_A} = k_1[A][E_A]-(k_2+k_3)[AE_A]$$

Equations at equilibrium

The equations described above are called mass action equations and their rate constants are not widely documented. What are actually documented are constants that describe a reaction scheme under the assumption that all enzyme complexes are at steady state (the rate of change is zero). In doing this you are able to remove enzyme and enzyme complex concentrations entirely from your equations. This is ideal as these concentrations are hard to measure in the lab.

In most cases these equations have been derived and documented especially when you note the equations are modular (you can chop the reaction scheme up into a series of smaller reaction schemes). If you have an abnormal reaction pathway you may need to derive it for yourself.

Concentration time graphs

The experimental data available to us was concentration time graphs, ODEs are perfect for creating these. We mainly used a matlab solver which is a predictor-corrector method. We also have a n^th order ODE solver written in c++ and that uses RK4.

It will now be explained how knowing the rate of change of your concentrations can lead to concentration over time in practice via such methods. The Forward Euler method will be used to illustrate.

You must know your starting concentrations, these are your “initial conditions” and represent e.g. your starting concentration of glucose/alcohol (mM). The differential equation is used to find the rate of change at this point (here on gradient). You take a small step in time along this gradient. Then you recalculate the gradient and take another small step and so on. See below figure for a visual representation.

Figure representing the euler method the dots represent recalculation points

A question we tried to answer is what reaction scheme is our cell free mechanism actually following. What follows are the possibilities we investigated, the reaction scheme, the equations for said reaction scheme and what they say physically about the reaction scheme. How to then predict which model/models are closest to the truth is explained in the network mechanism analysis.

Before we proceed it’s important to know that our schemes are broken up in two steps. Step 1 alcohol/glucose as a substrate with H₂O₂ as a product and step 2 H₂O₂ as a substrate with oxidised ABTS as a product. The final pairings we tested are shown in the chart below any pairing from each column is suitable e.g. Reversible Michaelis-Menten and Ping-Pong. The reason for these distinctions apart from saying Irreversible Michaelis-Menten is too far removed from reality is that they make no sense as the other step this will soon become clear.

Type of reaction	GOx reaction	HRP reaction
Irreversible Michaelis-Menten	n	n
Reversible Michaelis-Menten	y	n
Bi-Uni Michaelis-Menten	n	y
Uni-Bi Michaelis-Menten	y	n
Ping-Pong	n	y

For all the equations and variable definitions follow this link, Note for Tau = [P]/[S] for multiple substrates/products it is multiplicative. V represents the rate of change of the product.

Irreversible Michaelis-Menten. This model has no real physical relevance to our experiment so it's not discussed here. While useful as an initial test equation the irreversible nature easily hides properties of the system and this is highly undesirable

Reversible Michaelis-Menten:

This reaction scheme is saying alcohol/glucose and AOx/GOx can combine to make H₂O₂ and that this reaction can happen in reverse and there will be some sort of equilibrium.

The mass action scheme is:

As you can see the two schemes are very similar, we know that Ethanal is actually the second product in this reaction. Our question is, does it make a difference to the kinetics? This question can be assessed with the use of ensemble outputs and comparison to results using the Uni-Bi Michaelis-Menten rate law.

Uni-Bi Michaelis-Menten:

This reaction is saying that alcohol/glucose and AOx can combine to make H₂O₂ and Ethanal and that this can occur on reverse.

The mass action scheme is:

This mechanism does include Ethanal. The order the substrates break off does matter, as you can imagine if Ethanal does break off first you will end up with a large amount of Ethanal forcing the back reaction, this would be alleviated if H₂O₂ could break off first as then H₂O₂ would be more available for the next use. In literature there are two equations for this, however the ordered case is the random order case in the limit some of its rates go to zero. We assumed that if this were the case the sampling of our system would account for this, sampling many tiny parameters, rather than us making that judgement ourselves.

Bi-Uni Michaelis-Menten:

This reaction is saying H₂O₂ and ABTS and HRP can combine to make oxidised ABTS with HRP acting as a catalyst.

The mass action scheme is:

This mechanism takes into account that there are two substrates the same implentation was used as in Uni-Bi Michaelis-Menten. Physically it is saying there is some mechanism by which HRP acts directly as a catalyst in one step. Ping pong proposes something different. We can investigate which one is correct with our ensemble outputs.

PingPong:

This reaction is saying H₂O₂ and HRP can combine oxidising HRP and reducing H₂O₂, The resultant H₂O molecule takes no further place in this reaction but now the ABTS does. In this next step the reaction is saying ABTS and HRP* react reducing HRP and oxidising ABTS, The result is oxidised ABTS and HRP being reformed as a catalyst.

The mass action scheme is:

There is another alternative for a ping pong scheme and that is HRP and ABTS combine first oxidising the ABTS and reducing the HRP. The H₂O₂ would come along later and reform the HRP catalyst. This wasn’t considered as you would expect to see a jump in ABTS_oxidised at the start of your reaction and this wasn’t observed.

Practical considerations

Solving ODEs code

As mentioned in the theory section matlab has its own solver, ode45 was the one we used googling ode 45 will take you to MATLAB's documentation on it you can navigate to all the solvers it has. In these there is a concept called stiffness which you should understand to make an informed decision. However ode45 will work for enzyme kinetics with sensible initial conditions so it's not discussed here.

The c++ code uses the RK4 method. Its a nth order solver. All you have to do to use is enter your initial conditions into the first line of the code and in the diffall finction you must where there is a list starting with i==1 followed by an equation, have as many i==n statements as initial conditions with the equation written below corresponding to the nth initial conditions differential equation.

Note that all the above assumes your using first order differential equations only. However if you have a nth order differential equation it can be written as n first order equations. This mathematical step will need to be done first.

Discussion of why we chose our 4 models

When modelling the cell free system of the alcopatch there are a bank of potential rate equations, each with a differing underlying mechanism, for each step of the two steps.

GOx reaction	HRP reaction
Irreversible Michaelis-Menten	Irreversible Michaelis-Menten
Reversible Michaelis-Menten	Reversible Michaelis-Menten
Random Uni-Bi Michaelis-Menten	Random Bi-Uni Michaelis-Menten

This many potential equations mean there are 9 different combinations of rate laws More generally the this can be expressed as : $$n_{permutation}=\prod_{i=1}^{n_{steps}}n_{equations}$$

This can very quickly add to to huge numbers, a 5 step reaction pathway with 5 potential equations per step has 3125 possible combinations. The first step when building the reaction pathway is to eliminate any infeasible reaction schemes, just removing one potential reaction per step reduces this number by 2101! A short amount of critical consideration to remove equations is much better spent than trying to brute force a larger number of equations.

We chose four as they encapsulate the two questions we are asking, does Ethanal/GDL make a difference? And does the concentration of ABTS have aneffect? These four are irreversible Michaelis-Menten both steps, reversible Michaelis-Menten both steps. Reversible Michaelis-Menten followed by Random Bi-uni and Random uni-bi followed by Random Bi-uni. Ping pong was not used due to time constraints.

Extension to ensemble modelling

Extending to ensemble modelling is beautifully simple, You simply loop through the process multiple times using different rate constants generated from your pdf’s and save all the resulting concentration vs time arrays.

Constraints

Improve data motivation, What we used

In some situations there are constraints (representing physical limitations) that can be applied to the parameters to further improve your data set. One such example is the haldane relationship which relates kinetic parameters to the equilibrium constant.

For our system we noticed the Kcat/Km ratio was reported so this was used.

As such only parameters consistent with this ratio were used.

In doing this kind of step you not only save computer time but also can provide a tighter prediction on what you would expect to observe if the model is accurate.

Below you will find a figure showing the difference to the simulation using a filtered data set and a unfiltered data set.

The above figures illustrate the difference in simulation results when applying constraints to your parameters. Left is a concentration graph with no constraints applied, Right is a concentration graph with constraints applied. As you can see the constrained simulation produces less results that take much longer to reach steady state, this was consistent with experimental data.

<

Illustration of concept

Below some imaginary constraints are given by the 2 blue curves. Any point which falls between the 2 curves would be accepted (in green) and any point which lies outside these bounds is rejected (in red). This is repeated until enough suitable points are generated.

How code works and Git

Link to github homepage

The code works not exactly like you would expect from the theory explanation. Specifically the constraint is not applied to the original pdfs to create a new pdf consistent with the constraint, to do this you need to look into multivariate distributions. We actually took the much easier option of while running the simulation sampling from both Kcat and Km distributions (kcat/km is the constraint). Only keeping values that are within 0.5 standard deviations of the kcat/km value where standard deviation was calculated from the literature ratio.

The downside of our simpler method is that we waste computer time. As we only had one constraint this wasn’t much of a problem, however if you have many constraints you will have to look into multivariate techniques, or overcome the slowness with high performance computing (supercomputers).