This is my first post ever on this forum. I am a medical doctor working within the field of clinical pharmacology, and not a mathematician.
Bi- or triexponential decay functions are often used to model how drug concentrations change with time. In pharmacology, the bi-exponential decay function is typically formualted as follows:
$C(t)=A\cdot e^{-k_1\cdot t}+B\cdot e^{-k_2\cdot t}$
And the tri-exponential decay function:
$C(t)=A\cdot e^{-k_1\cdot t}+B\cdot e^{-k_2\cdot t}+D\cdot e^{-k_3\cdot t}$
$t$ is the time in hours. $C(t)$ is the drug concentration at time $t$. $k_1$, $k_2$ and $k_3$ are rate constants. $A$, $B$ and $D$ are coefficients whose sum is the initial drug concentration at time zero. All parameters are positive, real numbers.
I am working on a digital tool that requires the following question to be answered: how many hours must pass until the drug concentration is $x$ (where $x$ is some arbitrary number)?
Is it possible to find the inverse function of the bi- and tri-exponential functions above? If not, is there any other way for to easily solve the problem without using graphical or iterative methods?
Any help is much appreciated. Thank you.
Best Answer
In the most general case, you cannot inverse the sum of exponential functions and you need some numerical method (Newton being the simplest.
But, instead of trying to solve for $t$ $$C(t)=\sum_{i=1}^n a_i\, e^{-k_i t}$$ solving $$\log[C(t)]=\log\Big[\sum_{i=1}^n a_i\, e^{-k_i t}\Big]$$ would be much better since the rhs is much closer to linearity.
So, let us say that we look for the zeo of function $$f(t)=\log\Big[\sum_{i=1}^n a_i\, e^{-k_i t}\Big]-k$$ $k$ being given.
We have $$f(0)=\log\Big[\sum_{i=1}^n a_i\Big]-k >0\qquad \text{and}\qquad f''(0) >0$$ So, by Darboux theorem, starting iterations at $t_0=0$ ensures that we shall never face an overshoot of the solution.
But we can try to make an estimate trying to write $$\sum_{i=1}^n a_i\, e^{-k_i t}\sim A\,e^{-K t}=C(0)\,e^{-K t}$$ and define $$\frac {\sum_{i=1}^n a_i }K=\sum_{i=1}^n \frac {a_i}{k_i}\implies K=\frac{\sum_{i=1}^n a_i } {\sum_{i=1}^n \frac {a_i}{k_i} }$$ to match the areas under the curves from $t=0$ to $\infty$.
Trying withe random numbers $$a_1=29.8488\qquad a_2=50.0334\qquad a_3=21.9958$$ $$k_1=0.109865\qquad k_2=0.0760897\qquad k_3=0.0650516$$ gives $A=101.878 $ and $K=0.0803852$.
By curiosity, plot the two curves (they almost overlap).
Suppose that the target is $C(t)=45.678$; this gives as an estimate $t=9.97895$.
Now, Newton iterates $$\left( \begin{array}{cc} m & t_m \\ 0 & 9.97895 \\ 1 & 9.76084 \\ 2 & 9.76091 \end{array} \right)$$
Edit (for the biexponential decay function)
Since it is better to work on a logarithmic scale, consider that we want to approximate
$$\log\Big[ a_1\, e^{-k_1 t}+ a_2\, e^{-k_2 t}\Big]\qquad \text{by} \qquad \log\big[ a_1+ a_2\big]- K t$$ and consider two points $t_i=\frac{\log(2)}{k_i}$ corresponding to the half-times.
Minimzing the sum of the squared errors leads to $$(t_1^2+t_2^2) K=(t_1+t_2) \log(a_1+a_2)-$$ $$\Big[t_1 \log \left(a_1 e^{-k_1 t_1}+a_2 e^{-k_2 t_1}\right)+t_2 \log \left(a_1 e^{-k_1 t_2}+a_2 e^{-k_2 t_2}\right)\Big]$$
Trying with $a_1=30$, $k_1=0.11$, $a_2=50$, $k_2=0.07$, this would give $K\sim 0.0834$. Trying for a few values of $C(t)$
$$\left( \begin{array}{ccc} C(t) & \text{estimate} & \text{solution} \\ 70 & 1.60118 & 1.57641 \\ 60 & 3.44960 & 3.40984 \\ 50 & 5.63582 & 5.59716 \\ 40 & 8.31154 & 8.30198 \\ 30 & 11.7611 & 11.8338 \\ 20 & 16.6231 & 16.8954 \\ 10 & 24.9346 & 25.7633 \end{array} \right)$$
Using this estimate, Newton method would converge very fast. For example, using the last (and worst) point, we would have $$\left( \begin{array}{cc} n & t_n \\ 0 & 24.9346 \\ 1 & 25.7623 \\ 2 & 25.7633 \end{array} \right)$$