How Simple Numerical Integration Can Make Your Life Easier in Equation Modeling Tasks
Simulations and numerical modeling
In the Mathematics of natural sciences, engineering, and economics, the art of numerical integration serves as a powerful tool for describing the behavior of dynamic systems over time. The technique has a huge number of concrete applications to many problems across all these fields, especially when dealing with differential equations.
Technically, the main challenge arises when differential equations mix up one or more variables together with their derivatives in a way that analytical integration ends up entangling the different occurrences of the variables in such ways that they can't be put together and solved as a function of the other variables. In such cases, numerical integration, rather than regular integration, is at the rescue.
Although there are whole studies and methods to achieve numerical integration of differential equations efficiently, such as the Runge-Kutta methods, we can go to the basics and tackle the problem directly in a very simple "Euler-like" way, that is by using the instant values of the derivative and integrating small changes over time. This is exactly the case of the example application I will show you here: numerical integration of the Michaelis-Menten equation that describes enzyme catalysis. All with full maths, a working web app, and all its simple JavaScript code.
This article / tutorial is similar in style to this other one that I shared recently on how seemingly difficult procedures are actually rather simple and extremely helpful:
The power and simplicity of propagating errors with Monte Carlo simulations
Into differential equations and their simple, direct numerical integration
Differential equations often emerge when modeling dynamic systems, where variables change with respect to one another and over time. Such equations frequently involve derivatives of time, representing rates of change. However, their inherent complexity can make it difficult to extract intuitive insights or directly calculate specific outcomes, such as predicting when a system will reach a particular state.
As they are, the main issue with differential equations is their failure to provide a straightforward relationship between time and system variables. They describe how variables change continuously, but they don't readily yield insights into specific points in time when particular conditions or concentrations might be reached. This limitation can hinder our ability to make precise predictions and understand the dynamics of the systems under consideration.
The obvious way to sort out this problem is integrating the terms that contain derivatives. But often in the sciences, when you integrate a variable that already appears in the equation it ends up separated in different terms that can't be grouped together to solve for the variable. You will see this in the example of the Michaelis-Menten equation for enzyme catalysis that I discuss here.
An alternative is to use some numerical methods, among which the so-called Runge-Kutta methods are especially designed to solve these kinds of problems. But there's a much simpler numerical method, that works directly and explicitly with the very definition of derivatives: small change in one variable due to the small change in another, such as a position over time, a concentration over time, etc.
By discretizing time and approximating the continuous changes in variables, very simple numerical integration methods allow us to calculate how a system evolves over time in small, manageable steps. These methods enable us to approximate the behavior of dynamic systems at discrete time points, making it possible to predict when specific conditions will be met or to analyze the system's behavior at distinct moments.
Here I will show you how numerical integration can very easily make a complex equation like the Michaelis-Menten equation tractable and useful in a direct form, as a fundamental model in enzymology. This equation, which describes the kinetics of enzyme-catalyzed reactions, provides an ideal case study to showcase the elegance and versatility of numerical integration, which essentially breaks down complex dynamics into manageable steps. Moreover, the analytically integrated form of this equation presents the challenge that its most relevant variables can't be solved for each other, thus making the numerical integration method even more relevant.
Setting the example problem
The Michaelis-Menten equation is a fundamental model in enzymology that describes the rate of enzymatic reactions involving a substrate (S) and an enzyme (E) to produce a product. This equation plays a pivotal role in understanding the kinetics of enzyme-catalyzed reactions. However, obtaining an explicit solution for the substrate concentration as a function of time is impossible because there's no way of moving its terms around that will get you the concentration of substrate expressed as an analytical function of time.
But, as I will show you, we can numerically integrate the equation to simulate substrate concentration over time given some starting substrate and enzyme concentrations as well as a set of kinetic parameters. I will also put all this into work in the form of a web-based tool you can use to run your own simulations of substrate concentration over time very easily, just in your web browser, and whose code you can consult and copy for your ad hoc modifications.
The Michaelis-Menten Equation
The Michaelis-Menten equation is expressed as:
where:
V is the reaction velocity in M/s (molar per second) at a given concentration of enzyme and substrate, [E] is the enzyme concentration, typically talked in nano to micromolar but expressed here in molar units for consistency, [S] is the substrate concentration also in M, Km is the Michaelis constant describing the substrate binding part of the catalytic cycle (here also in M) and _k_cat is the turnover number or catalytic constant (in 1/s).
The equation describes how the reaction velocity V, which is the derivative of the substrate concentration [S] according to time t, depends on the concentrations of enzyme and substrate given a set of catalytic parameters _k_cat and Km.
Implicit variables, here substrate concentration
While the Michaelis-Menten equation provides a valuable insight into enzyme kinetics, it cannot be solved for the substrate concentration as an explicit function of time (which appears only upon integration of V as the derivative of [S] with respect to t). At the core, the limitation arises from the nonlinear and differential nature of the equation.
However, we can overcome this limitation through numerical methods, in particular by numerically integrating the reaction velocity which can be though as deltaConcentration/deltaTime when deltaTime is very small, that is d[S]/dt.
Numerical integration is for many problems a very powerful mathematical technique that approximates the behavior of a dynamic system over time by breaking it into small time steps during which something happens and variables change.
Let's see how it works here:
- To begin, we initialize the system with known values: enzyme concentration ([E]), initial substrate concentration ([S]0), catalytic constant (_k_cat), and Michaelis constant (Km). Moreover, we define a time step (dt) for the integration of the rate of change of substrate concentration, d[S]/dt = V. After a time step we will get the concentration [S] at t = tprevious + dt, and the total time length through which we track substrate concentration is given by a number of iterations, that we also must set.
- At each iteration, we calculate the instantaneous rate of change of [S] due to enzymatic reaction, which is represented by the reaction velocity V computed with the Michaelis-Menten equation just as you saw it above, plugging in the concentration of substrate that is left to be consumed, [S].
- This instantaneous velocity describes how fast the reaction is proceeding at a given moment. We can think that at that point of the reaction, in an infinitesimal amount of time then a proportional amount of substrate will be consumed. Since we can't have an infinitesimal number in the computer, we approximate this to some very small number (that we might need to tune, comments about this later). We then multiply the velocity calculated in step 2 by our deltat, or timestep, to determine the change in substrate concentration that would occur during this short time interval. This represents the amount of substrate that gets consumed (or produced) within this tiny time step.
- We then update the substrate concentration by subtracting its change as calculated in the previous step (V x dt) from the current concentration [S]. This step reflects the actual progression of the reaction over the short time interval, and gives as the new [S] at the current time.
- We repeat this process for a specified number of iterations, each time recalculating the velocity based on the updated substrate concentration.
- Once a series of [S] was calculated at increasing times t, we get the profile we were after.
Achieving accurate integration
The choice of the time step (dt) is crucial in numerical integration. A smaller time step allows us to capture finer details of the reaction dynamics but requires more computational resources. Conversely, a larger time step speeds up calculations but may lead to less accurate results or even break the whole thing -in this example by giving unrealistically large changes in substrate concentration thus introducing substantial biases or even producing such large changes that they are larger than the current substrate concentration itself, thus resulting in negative substrate concentrations.
Therefore, in practice, users should fine-tune the time step to strike a balance between computational efficiency and accuracy. The best in my experience is to try to find a small number that gives a curve with a reasonable shape resembling what's expected (see examples below) and then slowly reduce or increase the dt while checking if the shape of the curve is affected or not. If it's not, then your numerical integration is going well.
Example simulations with an online web app
This iterative process of calculating instantaneous velocity, updating substrate concentration, and repeating for many small time steps allows us to simulate the dynamic behavior of the reaction, even when an explicit mathematical solution is unattainable. And it's super easy to make a program out of this, just like in the web app that you can try here:
Here's an example simulation:
Such a tool allows users to explore the dynamic behavior of enzyme-catalyzed reactions very easily. In particular, understanding how substrate profiles change with different parameters helps to design experiments; for example how much enzyme to add to a reaction, how long one can expect to wait till a given fraction of substrate is consumed, etc. All things relevant in every-day studies of enzyme properties at laboratories doing from fundamental chemical biology to pharma companies. And all of them assisted by maths and computer code!
Some code—although, really, it's very simple
I find it beautiful that such a complex mathematical problem can be solved with a strategy like this one, that uses brute force but can be coded so simply.
While I invite you to see the full code of my web app with Ctrl-U in most web browsers, here's its core:
var lastS = S0 //Set to the starting concentration of substrate
var time = 0
for (var i = 0; i <= iterations; i++) {
var Vinst = kcat * E * lastS / (Km + lastS) //Instantaneous velocity at the current substrate concentration
var deltaS = Vinst * dt //Integration by multiplying the velocity by a sufficiently small timestep
lastS = lastS - deltaS //Update substrate concentration
time = time + dt //Update total time
data.addRow([time, lastS, S0 - lastS]); //Data is then plotted using Google's chart library in my web app
}Comparison to more complex methods for solving differential equations numerically
More complex methods specifically designed to solve these and other related kinds of problems exist, the main ones consisting in the Runge-Kutta series of methods. While I could cover Runge-Kutta methods here in the future, for the moment let me share with you this dedicated Wikipedia article and a subsequent comparison to the method I proposed here:
The method I presented above does not explicitly implement a specific Runge-Kutta method, but it does resemble a numerical integration scheme for solving ordinary differential equations. The numerical method we saw is much simpler and has the advantage that it conceptually relates very directly to the mere definition of derivative, here velocity as a "small change in a small amount of time", and using basic arithmetics.
One important disadvantage of the simple method we used is that it is sensitive to time step, not only if this parameter is too long but also if the velocities involved are too high. In the web app, this limitation is somehow compensated by the option to change timestep, with which you can repeat the simulation to easily check if the simulation is robust.
Compared to Runge-Kutta methods, which are designed to converge to the true solution as the time step approaches zero, the method used here may require such small time steps to achieve acceptable accuracy, that it might make the calculation computationally too expensive (because smaller time steps require more iterations to achieve long simulations).
For application of Runge-Kutta and even other different methods to the problem of the Michaelis-Menten equation, check out these papers:
Substrate and enzyme concentration dependence of the Henri-Michaelis-Menten model probed by…
Solution of the Michaelis-Menten equation using the decomposition method – PubMed
Takeaways
I hope I've convinced you here that numerical integration is a simple yet powerful mathematical technique that enables to solve problems like this one where it's impossible to get an explicit solution for the integrated variable. As you saw, the procedure breaks down the complex dynamics into manageable steps with an Euler-like idea rooted in the very concept of derivatives, thus being very simple from the viewpoint of the underlying maths and of the Programming involved.
www.lucianoabriata.com I write and photoshoot about everything that lies in my broad sphere of interests: nature, science, technology, programming, etc. Subscribe to get my new stories by email. To consult about small jobs check my services page here. You can contact me here.
