I wonder if anyone could please help me with this problem.
I have the following system of 2 nonlinear differential equations:
$$\frac {dM} {dt} = k_1 \cdot M \cdot (k_2 \cdot D - k_3) - k_4 \cdot M$$ $$\frac {dD} {dt} = - k_1 \cdot M \cdot (k_2 \cdot D - k_3) - (k_4 + k_5) \cdot D$$
I know the value of $M$ and $D$ at time $t=0$; all constants $k_i$ are real and positive.
To give some context, the system is a simplified description of gastrointestinal (GI) disposition of a drug. $M$ is the amount of drug present as a solid, and $D$ is the total amount of dissolved drug. The initial term in the rhs of each equation (containing $k_1, k_2, k_3$) is the Nernst-Brunner dissolution-precipitation term; as you can see, if there were no other terms, an equilibrium would be reached between solid and dissolved drug. The other terms are related to excretion from the GI tract ($k_4$), which happens to both the solid and dissolved drug equally, and the sum of absorption and GI decomposition ($k_5$), which only operates on $D$.
Ideally I would like to find explicit expressions for $M(t)$ and $D(t)$, but as far as I know, the usual techniques I would use (like Laplace) don't work in this case. I also tried taking the ratio between the two equations, or substituting $R=M+D$, because the sum of the two equations doesn't contain any nonlinear terms, but I don't seem to be able to find a solution.
Q1: do you think there is an analytical solution at all, and if so, could you please advise how to obtain it?
Lacking that, I would be OK with a mass balance.
In the physical system I described, $M$ and $D$ go to $0$ when $t$ goes to infinite. Therefore, if I multiply both equations by $dt$ and integrate, I know the lhs, and I can leave the integrals of $M$ and $D$ implicit:
$$0-M(0) = k_1 \cdot k_2 \cdot \int_0^\infty M \cdot D \cdot dt - (k_1 \cdot k_3 + k_4) \cdot \int_0^\infty M \cdot dt$$ $$0-D(0) = - k_1 \cdot k_2 \cdot \int_0^\infty M \cdot D \cdot dt + k_1 \cdot k_3 \cdot \int_0^\infty M \cdot dt - (k_4 + k_5) \cdot \int_0^\infty D \cdot dt$$
In a linear system, at this point I would solve for the integrals of $M$ and $D$, which would then allow me to calculate the total amount of drug excreted, for instance:
$$D_{excreted} = k_4 \cdot \int_0^\infty (M + D) \cdot dt$$
and the sum of drug absorbed and decomposed:
$$D_{abs+dec} = k_5 \cdot \int_0^\infty D \cdot dt$$
which also respects the overall mass balance:
$$D_{excreted}+D_{abs+dec} = ... = M(0)+D(0)$$
But I can't do any of that, because I have the integral of $M \cdot D$ in both equations and I don't know how to get rid of it.
Q2: can you suggest any way to get at least the mass balance, i.e. an explicit expression of $D_{abs} ?$
Thanks
EDIT (after posts by player100)
By eliminating the unwanted integral of $M \cdot D$ from the system of the two mass balance equations, and solving for the integral of $D$, I get:
$$\int_0^\infty D \cdot dt = \frac {M(0)+D(0)-\int_0^\infty M \cdot dt} {k_4+k_5}$$
I tried to figure out if the integral of M could be somehow rewritten as a function of the integral of D, given player100's suggestion, but I'm stuck at:
$$M(t)=M(0) \cdot \frac {e^{k_1 \cdot k_2 \cdot \int_0^t D({\tau}) \cdot d{\tau}}} {e^{(k_1 \cdot k_3 + k_4) \cdot t}}$$
So my next move was to give up on the original Nernst-Brunner model, and use a simpler model with linear terms only:
$$\frac {dM} {dt} = k_P \cdot D - k_D \cdot M - k_4 \cdot M$$ $$\frac {dD} {dt} = - k_P \cdot D + k_D \cdot M - (k_4 + k_5) \cdot D$$
which of course can be solved both leaving the integrals implicit and also analytically for $M(t)$ and $D(t)$ (e.g. by Laplace).
As for the values of $k_P$ and $k_D$, I still need to find the best strategy to calculate them from the Nernst-Brunner parameters $k_1, k_2, k_3$.
The current approach I am trying is to consider the systems with no other terms but dissolution and precipitation, i.e.:
$$\frac {dM} {dt} = -\frac {dD} {dt} = k_1 \cdot M \cdot (k_2 \cdot D - k_3)$$
for the Nernst-Brunner (NB) case, and:
$$\frac {dM} {dt} = -\frac {dD} {dt} = k_P \cdot D - k_D \cdot M$$
for the simplified case.
They can both be solved analytically thanks to the fact that:
$$M(0)+D(0)=M(t)+D(t)$$
and it can be shown that $M$ and $D$ reach an 'equilibrium', i.e. a constant value, when t tends to infinite.
My aim is to choose $k_P$ and $k_D$ so that $D$ for the NB system is 'as close as possible' to $D$ for the simplified system.
To achieve that, one important condition is that the equilibrium values of $M$ (and $D$) must be the same in the simplified system and in the NB system. This gives me one equation in $k_P$ and $k_D$.
The other equation is what I'm still trying to figure out. At first I thought equating the rates at time $0$ would be OK, but the NB system is a 'self-catalysed' one, i.e. if I start with very little $M$ and the parameters are such that much more $M$ must be present at equilibrium, in the NB system the rate of formation of $M$ initially increases as more $M$ is formed, reaches a maximum and then decreases; instead, in the simplified system it always decreases, because $M$ becomes larger and $D$ becomes smaller. The opposite happens if there is much $M$ at the start and it must dissolve: there the rate is maximal at the beginning, but as $M$ decreases, it decreases much more rapidly than by the decrease of $D$ alone; the simplified system can't account for that. So, equating the initial rate means that the simplified system is too slow in forming solid when there is too little of it, and too fast in dissolving it when there is too much of it, compared to the more 'exact' NB case.
I already tried integrating the squared difference between the NB $D(t)$ and the simplified $D(t)$, but the integral is not analytic.
So I think I need to find a simpler strategy, for instance imposing that at the time when the value of $D(t)$ is 90% of its equilibrium value in the NB case, the same is true for the simplified case. That should give me the missing equation in $k_P$ and $k_D$.
Of course any suggestions regarding this new problem are welcome!
EDIT2
I found a way around the problem.
The integral of the squared difference between the NB $D(t)$ and the simplified $D(t)$, is indeed not analytic as I said, but the integral of their difference (not squared) is analytic, and as both curves have the same start and end value and are always $>=0$, there is no problem.
By setting the integral to $0$ I got my second equation and could solve the system for $k_P$ and $k_D$.
Numerical tests show that the match between the two models is not always great, but it's good enough for our purposes, and I don't think we could do much better with a linear model.
$\frac {d}{dt}ln(M) = k_1k_2D-k_1k_3$
$\frac {M}{M_0}= exp\left( \int_0^t{k_1k_2D-k_1k_3 \ d\tau} \right)$
$\begin{align} \frac{dD}{dt} &= -k_1M_0exp\left( \int_0^t{k_1k_2D-k_1k_3 \ d\tau} \right)(k_2D-k_3)-(k_4+k_5)D \\\ &= -k_1M_0exp\left( k_1k_2\int_0^t{D \ d\tau} \right)(k_2D-k_3)exp\left(-k_1k_3t\right)-(k_4+k_5)D \end{align}$
$e^{(k_4+k_5)t}\left(\frac{dD}{dt}+(k_4+k_5)D\right)=\frac{d}{dt}\left(e^{(k_4+k_5)t}D\right)$
$\frac{d}{dt}\left(e^{(k_4+k_5)t}D\right) = -k_1M_0exp\left( k_1k_2\int_0^t{D \ d\tau} \right)(k_2D-k_3)exp\left((k_4+k_5-k_1k_3)t\right)$
Let $E=e^{(k_4+k_5)t}D$, then
$\frac{d}{dt}E = -k_1M_0exp\left( k_1k_2\int_0^t{Ee^{-(k_4+k_5)\tau} \ d\tau} \right)(k_2Ee^{-(k_4+k_5)t}-k_3)exp\left((k_4+k_5-k_1k_3)t\right)$
This is amenable to numerical methods.