landercy.be

The Lambert function

Written by

jlandercy

in

,

The Lambert function [1783] is a less known function which has a lot of nice applications in sciences. This post presents its application in Biochemistry and Physical Chemistry in real situations.

The Lambert function WW is defined as follow:

z=w⋅exp⁡(w)⇔w=W(z)z = w\cdot\exp(w) \Leftrightarrow w= W(z)

It is a multi-valued function with two real branches W0W_0 and W−1W_{-1}. The figure below shows the real valued direct function and its inverses functions (real branches of Lambert function).

The Lambert function is also referenced as Product Log in some articles and softwares.

Michaelis-Menten

It can be shown [1997] that the Michaelis-Menten equation can lead to an analytical solution using the Lambert function. Given the Michaelis-Menten equation:

dSdt=−vmax⁡⋅SKM+S\frac{\mathrm{d}S}{\mathrm{d}t} = -\frac{v_{\max} \cdot S}{K_M + S}

It is a separable ordinary differential equation (ODE), we can solve it explicitly for tt:

∫S0S(KMvmax⁡S+1vmax⁡)dS=−∫0tdt\int\limits_{S_0}^{S}\left(\frac{K_M}{v_{\max}S} + \frac{1}{v_{\max}}\right)\mathrm{d}S = -\int\limits_{0}^{t}\mathrm{d}t

It goes as follow:

−t=Kvmax⁡ln⁡|SS0|+1vmax⁡(S−S0)-t = \frac{K}{v_{\max}}\ln\left|\frac{S}{S_0}\right| + \frac{1}{v_{\max}}(S-S_0)

This is where the Lambert function becomes useful, it allows us to isolate S(t)S(t). First we rewrite the equation such as the LHS has the form of the Lambert definition:

SKMexp⁡(SKM)=S0KMexp⁡(−vmax⁡t+S0KM)\frac{S}{K_M}\exp\left(\frac{S}{K_M}\right) = \frac{S_0}{K_M}\exp\left(\frac{-v_{\max}t + S_0}{K_M}\right)

Then we simply apply the definition of the Lambert function:

S(t)=KMW[S0KMexp⁡(−vmax⁡t+S0KM)]S(t) = K_MW\left[\frac{S_0}{K_M}\exp\left(\frac{-v_{\max}t + S_0}{K_M}\right)\right]

Which is the analytical solution of the Michaelis-Menten equation.

Numerical solution

The Michaelis-Menten ODE can be solved numerically using Python and scipy for a given set of parameters and initial condition:

import numpy as np
from scipy import special, integrate, optimize
import matplotlib.pyplot as plt

K_M = 1e-3
v_max = 1e-4

def system(t, S, K_M, v_max):
    return np.array([
        - v_max * S[0] / (K_M + S[0])
    ])

t = np.linspace(0, 300, 200)
S0 = 1e-2

solution = integrate.solve_ivp(system, [t[0], t[-1]], [S0], t_eval=t, args=(K_M, v_max))
S = solution.y[0,:]

Now we can recover this parameters and initial condition by fitting the ODE solution using the analytical function found above:

def model(t, K_M, v_max, S_0):
    return np.real(K_M * special.lambertw((S_0 / K_M) * np.exp((- v_max * t + S_0) / K_M), 0))

popt, pcov = optimize.curve_fit(model, t, S)
Shat = model(t, *popt)

We can show that both model numerically agree:

Analytical solution

The Michaelis-Menten ODE can be solved analytically using Computer Algebra System (CAS) such as Python and sympy package:

import sympy as sp

K = sp.Symbol("K", real=True, positive=True)
v = sp.Symbol("v", real=True, positive=True)
S0 = sp.Symbol("S0", real=True, positive=True)
t = sp.Symbol("t", real=True)
S = sp.Function("S")

equation = sp.Eq(sp.Derivative(S(t), t), - v * S(t) / (K + S(t)))

sp.dsolve(equation, ics={S(0): S0})

Which returns the above found function:

S(t)=KW(eKlog⁡(S0)+S0−tvKK)\displaystyle S{\left(t \right)} = K W\left(\frac{e^{\frac{K \log{\left(S_{0} \right)} + S_{0} – t v}{K}}}{K}\right)

It can also be confirmed using Wolfram CAS.

Heterogeneous Catalysis

With the assumption of the Langmuir Adsorption model, for a competitive heterogeneous catalysis such as deshydrogenation of isopropanol in presence of Raney’s Nickel catalyst. A typical ODE model of the reaction coordinate ξ\xi is:

v=1Vdξdt=kb1(n0−ξ)b1(n0−ξ)+b2ξv = \frac{1}{V}\frac{\mathrm{d}\xi}{\mathrm{d}t} = k \frac{b_1 (n_0 – \xi)}{b_1(n_0-\xi)+b_2\xi}

Where VV and n0n_0 are constants of the experimental setup and k′=kVk’=kV and α=b2b1\alpha=\frac{b_2}{b_1} are reaction constants, respectively the molar kinetic constant and the adsorption affinity constant ratio of product and reactant.

This ODE is separable:

∫0ξ(1+αξn0−ξ)dξ=k′∫0tdt\int\limits_0^\xi\left(1 + \alpha\frac{\xi}{n_0 – \xi}\right)\mathrm{d}\xi = k’\int\limits_0^t\mathrm{d}t

Integration w.r.t. initial conditions leads to:

ξ(1−α)+αn0log⁡(n0−ξn0)=k′t\xi (1 – \alpha) + \alpha n_0 \log\left(\frac{n_0 – \xi}{n_0}\right) = k’t

Which is an explicit function t(ξ)t(\xi). For α=1\alpha=1 the equation is easily invertible. For α≠1\alpha \neq 1, it is also invertible by the mean of the Lambert function, it is then possible to get an analytical form of ξ(t)\xi(t) in both cases.

Recalling that ba=exp⁡(alog⁡(b))b^a = \exp(a\log(b)), it requires a bit of algebra to show that:

(n0−ξn0)(α1−α)exp⁡(n0−ξn0)=exp⁡(1−n0k′t1−α)\left(\frac{n_0 – \xi}{n_0}\right)^{\left(\frac{\alpha}{1 – \alpha}\right)}\exp\left(\frac{n_0 – \xi}{n_0}\right) = \exp\left(1 – \frac{n_0 k’ t}{1 – \alpha}\right)

Now we apply the following property (see next section):

waexp⁡(w)=z⇔w=aW[zaa]w^a\exp(w) = z \Leftrightarrow w = aW\left[\frac{\sqrt[a]{z}}{a}\right]

To get the expression of the reaction coordinate:

ξ(t)=n0(1−(α1−α)W[exp⁡(1−α−n0k′tα)α1−α])\xi(t) = n_0 \left( 1 – \left(\frac{\alpha}{1 – \alpha}\right)W\left[ \frac{\exp\left( \frac{1 – \alpha – n_0k’t}{\alpha} \right)}{\frac{\alpha}{1 -\alpha}} \right] \right)

The solution of the original ODE is then:

ξ(t)={n0(1−exp⁡(−k′n0t))if α=1n0(1−(α1−α)W[exp⁡(1−α−n0k′tα)α1−α])if α≠1\xi(t) = \begin{cases} n_0\left(1 – \exp\left(-\frac{k’}{n_0}t \right)\right) & \text{if } \alpha = 1 \\ n_0 \left( 1 – \left(\frac{\alpha}{1 – \alpha}\right)W\left[ \frac{\exp\left( \frac{1 – \alpha – n_0k’t}{\alpha} \right)}{\frac{\alpha}{1 -\alpha}} \right] \right) & \text{if } \alpha \neq 1 \end{cases}

Numerical Solution

We first solve the ODE system with scipy:

kp = 1e-3 * 200e-3
a = 70.
n0 = 2.5

def system(t, xi, kp, a, n0):
    return np.array([
        kp * (1. / (1. + a * xi[0] / (n0 - xi[0])))
    ])

t = np.linspace(0, 8000, 200)

solution = integrate.solve_ivp(system, [t[0], t[-1]], [0], t_eval=t, args=(kp, a, n0), method="LSODA")
xi = solution.y[0,:]

From the solution, we can derive parameters using Non Linear Least Square (NLLS):

def model(t, kp, a, n0):
    if np.abs(a - 1.) <= 1e-8:
        return n0 * (1. - np.exp(- kp * t / n0))
    else:
        return np.real(n0 * (1. - a / (1 - a) * special.lambertw(np.exp((1 - a - n0 * kp * t) / a) / (a / (1 - a)), 0)))

popt, pcov = optimize.curve_fit(model, t, xi, [1e-3, 1e2, 1.], bounds=[[0.] * 3, [np.inf] * 3])
xihat = model(t, *popt)

It required some fine tuning such as genuinely choosing the initial parameters and fixing bounds. Both solutions numerically agree:

Analytical solution

We can confirm our calculations using Python and sympy:

t = sp.Symbol("t", real=True, positive=True)
k = sp.Symbol("k", real=True, positive=True)
a = sp.Symbol("a", real=True, positive=True)
n0 = sp.Symbol("n0", real=True, positive=True)
x = sp.Function("x")

equation = sp.Eq(
    sp.Derivative(x(t), t),
    k * (1 / (1 + a * x(t) / (n0 - x(t))))
)
sp.dsolve(equation, ics={x(0): 0})
sol.rhs.limit(a, 1)

Which returns an equivalent solution:

x(t)=n0(aW((a−1)ean0log⁡(−n0)−an0−kt+n0an0an0)+a−1)a−1\displaystyle x{\left(t \right)} = \frac{n_{0} \left(a W\left(\frac{\left(a – 1\right) e^{\frac{a n_{0} \log{\left(- n_{0} \right)} – a n_{0} – k t + n_{0}}{a n_{0}}}}{a n_{0}}\right) + a – 1\right)}{a – 1}

And properly converges to the degenerated solution when taking the limit α→1\alpha \rightarrow 1:

n0(1−e−ktn0)\displaystyle n_{0} \left(1 – e^{- \frac{k t}{n_{0}}}\right)

Verification with Wolfram CAS gives same results.

Properties

Proof

Given the following equation:

z=waexp⁡(w)z = w^a\exp(w)

We take the aa-th root:

za=wexp⁡(wa)\sqrt[a]{z} = w\exp\left(\frac{w}{a}\right)

Dividing booth side by aa:

zaa=waexp⁡(wa)\frac{\sqrt[a]{z}}{a} = \frac{w}{a}\exp\left(\frac{w}{a}\right)

Finally applying the Lambert definition, we get:

W[zaa]=waW\left[\frac{\sqrt[a]{z}}{a}\right] = \frac{w}{a}

Which leads to the desired property:

w=aW[zaa]w = aW\left[\frac{\sqrt[a]{z}}{a}\right]

Origin

Another useful property is the following identify which holds for x>0x>0:

xW[exp⁡(1x)x]=1xW\left[\frac{\exp\left(\frac{1}{x}\right)}{x}\right] = 1

And is trivial to show just by applying the Lambert definition to 1x\frac{1}{x}. At the limit of zero, it diverges, one can show that:

limx→0−⁡xW[exp⁡(1x)x]=0\lim_{x \rightarrow 0^-} xW\left[\frac{\exp\left(\frac{1}{x}\right)}{x}\right] = 0

And:

limx→0+⁡xW[exp⁡(1x)x]=1\lim_{x \rightarrow 0^+} xW\left[\frac{\exp\left(\frac{1}{x}\right)}{x}\right] = 1

The following limits involves indeterminate forms of ∞∞\frac{\infty}{\infty}, they can be reduced by using the l’Hospital rule and the derivative formula of the Lambert function:

ddxW(x)=W(x)xW(x)+x\frac{\mathrm{d}}{\mathrm{d} x}W(x)=\frac{W(x)}{x W(x)+x}

This property ensures that the solution actually pass by the origin as required by the initial conditions.

Conclusions

We showed through applied chemistry examples that the Lambert function can provide analytical solution to common ODE systems.

The Lambert function comes into play after the integration part when it needs to get the inverse function of moderately complex expressions. The trick is always the same, performing algebra to express the product wexp⁡(w)w\exp(w) and then applying the Lambert function.

We also showed that Python equipped with scipy is an useful package to numerically solve the ODE or to implement regression of curves involving the Lambert function and that Python equipped with sympy package is helpful in solving the ODE and doing the necessary algebra to apply the Lambert function.

Just as the logarithm function is, the Lambert function is a useful addition to the analyst’s toolkit.

More posts